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ABSTRACT 

K  ■- 

This  thesis  presents  the  theoretical  and  computational 
underpinninas  of  a  novel  approach  to/ithe  determination  of'1  the 
acoustic  parameters  of  the  ocean  bottom  using  a  monochromatic 
source.  The  problem  is  shown  to  be  equivalent  to  that  of  the 
reconstruction  of  the  potential  in  a  Schrftdinger  equation 
from  the  knowledge  of  the  plane-wave  reflection  coefficient 
as  a  function  of  vertical  wavenumber,  for  all  real 

positive  kz*  First,  the  reflection  coefficient  is  shown  to  1 
decay  asymptotically  at  least  as  fast  as  Jl/kz2)  for  larqe  kz_ 
and  is  therefore  integrable.  The  Gelfand^Xevitan  inversion 
procedure  is  extended  to  include  the  case  of  basement 
velocity  hiqher  than  the  velocity  of  sound  in  water.  The 
peqlect  of  bound  states  is  shown  to  be  justified  in  both 
clayey  silt  and  silty  clay  at  the  220  Hz  frequency  of 
operation . 

Three  methods  for  the  numerical  solution,  of 1  the  integral 
eauation  are  investigated.  The  first  <one  is  an  "Improved 
Born  approximation"  wherein  the  solution  is  given  as  a  series 
expansion  the  fipst  term  of  which  is  the  Born  approximation 
while  the  second,' term  represents  a  substantial  and  yet  easy 
to  implement  improvement  over  Born. 

The  two  other  methods  are  based  on  a  discretization  of 
the  Gelfand-Levitan  inteqral  equation,  and  both  avoid  a 
matrix  inversion:  one  by  employinq  a  recursive  procedure, 
and  the  other  by  coupling  the  Gelfand-Levitan  equation  with  a 
partial  differential  equation.  Bounds  are  obtained  on  errors 
in  the  solution  due  either  to  discretization  or  to  data  inac¬ 
curacy.  These  methods  are  tested  on  synthetic  data  obtained 
from  known  qeoacoustic  models  of  the  ocean  bottom.  Results 
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are  found  to  be  very  accurate  particularly  at  the  top  of  the 
sediment  layer  with  resolution  of  less  than  the  wavelenqth  of 
the  acoustic  source  in  the  water.  Several  effects  are  inves- 
tiqated,  such  as  samplinq,  attenuation,  and  noise.  Also 
examined  is  the  qradual  restriction  of  the  reflection  coeffi¬ 
cient  to  a  finite  ranqe  of  vertical  wavenumbers  and  the  con¬ 
sequent  proqressive  deterioration  of  the  reconstruction. 

The  analysis  shows  how  to  reconstruct  velocity  profiles 
in  the  presence  of  density  variation  when  the  experiment  is 
conducted  at  two  frequencies. 

Our  results  provide  a  qood  understandi nq  of  the  issues 
involved  in  conductinq  a  monochromatic  deep  ocean  bottom 
experiment  and  constitute  a  promisinq  technique  for  process- 
ina  the  experimental  data  when  it  becomes  available. 
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CHAPTER  I 


BACKGROUND 


It  is  a  rare  person  whose  pulse  is 
not  stirred  by  the  dramatic  siqht  of 
the  restless  surface  of  the  sea. 

The  chaotic  sea  surface  is  a 
limitless  source  of  inspiration  to 
poet,  painter  and  musician  alike. 

But  what  lies  beneath  this  churninq 
surface?  How  can  we  probe  the 
depths  of  the  sea? 

C.  Clay  &  H.  Medwin 


1  .1  Background 

The  sea  floor  beqins  at  the  water-sediment  interface, 
overlies  the  sedimentary  layer,  and  beneath  it,  the  oceanic 
crust.  The  study  of  the  ocean  bottom  has  been,  until 
recently,  the  province  of  the  marine  qeoloqist  seekinq  to 
probe  the  oceanic  crust,  and  to  unravel  the  secrets  of  its 
structure  and  evolution.  The  marine  qeoloqist  has  now  been 
joined  by  the  underwater  acoustician  studying  the  transmis¬ 
sion  of  low  frequency  sound  throuqh  the  ocean;  It  has  become 
clear  in  liqht  of  underwater  sound  propagation  experiments 
carried  out  at  the  various  Oceanographic  Institutions ^ ^  that 
long  range  low  frequency  sound  transmission  is  affected  by 
the  nature  of  the  ocean  bottom,  and  hence,  that  acoustic  wave 
propagation  models  -  to  characterize  sonar  performance,  for 
instance  -  should  include  a  detailed  representation  of  the 


bottom,  the  qeoacoustic  model.  In  practice,  one  posits  a 
"reason able"  ocean  bottom  model,  and  then  one  proceeds  to 
?c‘.;e  the  propagation  problem  at  hand,  (the  "Direct" 
problem).  Here,  the  opposite  perspective  is  adopted:  Since 
the  ocean  bottom  affects  acoustic  wave  propagation,  would  it 
not  be  possible  to  learn  something  about  the  bottom  from  that 
interaction?  (The  "Inverse"  problem.) 

The  answer  to  this  question  is  being  sought  in  the 
context  of  an  original,  single  frequency,  deep  ocean  bottom 
interaction  experiment  desiqned  by  G.  Frisk  and  his 
colleaques  of  the  Woods  Hole  Oceanographic  Institution,  and 
performed  in  the  Hatteras  Abyssal  Plain^)»(3)  ancj  at  other 
locat ions ^ 4 ^ .  The  monochromatic  character  of  the  Frisk 
method  sets  it  apart  from  currently  used  techniques  usinq 
explosive  (or  impulse-like)  wide-band  sources.  The  Frisk 
experiment  started  off  as  a  heuristic  approach . 

This  thesis  presents  the  theoretical  basis  and  numerical 
analysis  of  the  monochromatic  experiment  based  on  an  exten¬ 
sion  of  the  Gelf and-Levi tan  theory  of  quantum  scattering.  We 
succeeded  in  applyinq  a  numerical  solution  to  the  exact 
inverse  method  which  distinguishes  this  solution  method  from 
the  currently  used  approximate  or  trial  and  error  inverse 
methods  . 
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1.2  The  Experiment 


The  aeonetry  of  the  experiment  is  best  explained  with 
reference  to  Fiq.  1.  In  brief,  a  drifting  vessel  tows  a 
220  Hz  pulsed  CW  source.  Hydrophones  moored  close  to  the 
bottom  record  the  resultinq  pressure,  both  amplitude  and 
phase  (via  coherent  quadrature  demodulation),  as  the  ship 
opens  ranqe. 

The  pulsinq  of  the  220  Hz  source,  turning  it  on  for 
4  sec  every  14  sec,  allows  for  steady  state  conditions  to  be 
attained  before  any  reflections  from  the  ocean  surface  can 
reach  the  receivers.  Note  that  a  14  second  duty  cycle 
siqnifies  that  the  acoustic  field  is  sampled  spatially  once 
every  half  wavelength. 

The  source  aperture  is  small  enouqh  compared  to  the 
wavelenqth  (7  m)  to  be  considered  an  omnidirectional  point 
source.  The  recorded  complex  pressure  is  therefore  the  field 
due  to  the  reflections  of  a  spherical  wave  off  the  bottom. 

The  information  is  translated  via  a  Hankel  transform,  into 
the  plane  wave  reflection  coefficient  at  a  single  frequency 
(220  Hz)  for  all  angles  of  incidence,  both  real  and 
complex^5).  The  critical  point  to  observe  is  that  for  a 
monochromatic  plane  wave  incident  on  a  flat  layered  bottom, 
the  reflection  coefficient  is  a  function  of  the  anqle  of 
incidence.  At  a  qiven  anqle  of  incidence,  the  maqnitude  and 
phase  of  the  reflection  coefficient  depend  on  the  acoustic 
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properties  of  the  bottom  sediments.  That  naturally  leads  to 
the  ocean  bottom  inverse  problem:  Can  the  sediment  acoustic 
parameters,  velocity  and  density,  be  reconstructed  from  the 
plane  wave  reflection  data  at  a  sinqle  frequency  for  all 
anqles  of  incidence?  The  thesis  aims  to  elucidate  that 
questions . 

1 . 3  Experimental  Data 

The  relationship  between  the  experimental  data  and  the 
plane  wave  reflection  coefficient  has  been  studied  numerical¬ 
ly  by  Frisk  et  al^^  and  Mook^).  Before  beqinninq  our 
analysis,  it  is  useful  here  to  take  a  quick  look  at  the 
underlyinq  theory  relatinq  the  reflections  of  spherical  waves 
from  a  point  source,  as  in  the  Frisk  experiment,  to  the  plane 
wave  reflection  as  in  our  model. 

The  measurements  yield  the  pressure  field  due  to  a  point 
source  above  the  bottom  half  space.  Because  of  the  cylindri¬ 
cal  symmetry  of  the  problem,  the  reflected  pressure  field  can 
be  written  as  a  superposition  of  plane  waves 


PR(r,z) 


(1.1) 

where  kr  is  the  horizontal  wavenumber,  and  r(kr)  is  the 
correspondinq  reflection  coefficient. 
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This  Sommerfeld  integral  is  a  Hankel  transform  which  can 


be  inverted  to  obtain  the  reflection  coefficient 


PR(r)J0(kr«Odr 


(1.2) 


Note  that  r (kr)  is  a  function  of  horizontal  wavenumber  while 
the  required  input  to  the  inverse  procedure  is  the  reflection 
coefficient  as  a  function  of  the  vertical  wavenumber  r(kz). 
The  two  are  of  course  related  by  the  dispersion  relation 


Given  r(kr)  for  real  kr,  one  can  readily  generate  r(kz)  for 
0  <  kz  <  kg.  It  is  more  difficult  to  obtain  r(kz)  for  the 
full  range  0  <  kz  <  ®  since  kg  <  kz  <  «  correspond  to  r (kr) 
for  imaginary  kr.  One  approach  has  been  suggested  by 
Stickler^  and  involves  the  use  of  a  theorem  by  Van  Winter 
to  generate  r(kz)  on  a  ray  in  the  complex  plane  given  its 
value  on  the  real  axis  segment  0  <  kz  <  kg.  The  effect  of 
limiting  r(kz)  to  real  angles  (0  <  kz  <  kg)  on  the  inversion 
for  the  unknown  potential  V(z)  will  be  discussed  in 
Chapter  VI. 
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1.4  The  Model 


A  number  of  assumptions  are  implicit  in  the  "exact" 
inverse  procedure  detailed  in  the  following  chapters. 

(a)  The  ocean  is  assumed  homogeneous,  and  acoustically 
transparent.  In  particular,  the  depth  variation  of  the 
velocity  of  sound  propagation  is  neglected.  That  is  a 
reasonable  assumption  at  the  great  depth  (5  km)  in  which  the 
experiment  is  conducted. 

(b)  The  ocean  bottom  is  assumed  to  have  no  horizontal 
structure,  the  velocity  variation  is  therefore  solely  a 
function  of  depth.  That  is  a  severe  restriction  imposed  by 
all  "exact"  inversion  formalisms  developed  to  date. 
Surprisingly,  horizontal  stratification  describes  adequately 
vast  areas  of  the  deep  ocean  floor  known  as  Abyssal  plains: 
These  are  widespread  in  the  Atlantic  and  Indian  Oceans  and  in 
the  marginal  seas. 

The  early  deep  ocean  bottom  interaction  experiments  were 
conducted  in  the  Hatteras  Abyssal  Plain.  This  nearly  level 
plain  lies  at  the  base  of  the  East  Coast  Continental  rise, 
and  is  1000  km  long  by  150-300  km  wide.  Its  thick  (>  1  km) 
sediments  were  formed  by  the  smooth  accumulation  of 
turbidites  over  the  rough  basement  resulting  in  one  of  the 
flattest  areas  on  earth  with  slopes  of  less  than  1  m/km. 
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(c)  The  ocean  bottom  has  been  traditionally  treated  as 
a  fluid  in  problems  involvinn  bottom  reflection.  In  the  deep 
ocean  bottom  interaction  experiments,  the  acoustic  source 
qenerates  compress iona 1  (P)  waves  in  the  water,  which  upon 
propagation  in  the  bottom,  a  vertically  hete roqeneous  medium, 
may  be  converted  to  shear  ( SV)  waves.  The  conversion  to 
shear  waves  will  occur  discretely  at  layer  interfaces  and 
continuously,  where  velocity  qradients  occur.  Fryer  had 
shown  in  one  of  his  papers^®)  that  coupling  for  continuously 
varying  elastic  parameters  is  negligible  at  frequencies  above 
20  Hz.  Vidmar  and  Foreman^)  estimated  that  gradient-induced 
couplinq  should  be  expected  in  marine  sediment  at  frequencies 
up  to  3  Hz.  Another  paper  by  Fryer ^ 10 ^  established  that  this 
coupling  is  extremely  small  above  1  Hz,  regardless  of 
sediment  thickness.  The  most  important  effect  of  coupling 
appears  to  be  the  conversion  of  shear  to  compressional  motion 
at  the  sediment  basement  interface.  Note,  that  although 
these  results  are  based  on  a  continuously  varying  structure 
(approximated  by  homogeneous  layers),  they  do  provide  for  the 
sharp  discontinuity  in  elastic  parameters  at  the  sediment 
basement  interface.  These  results  do  justify  the  neqlect  of 
shear  wave  effects  at  the  220  Hz  frequency  selected  for  the 
experiments  that  have  already  been  conducted,  and  at  the 
lower  frequencies  envisaged  by  the  Frisk  qroup  for  future 
experiments . 
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1.5 


Overview  of  the  Thesis 


In  Chapter  II,  we  begin  with  a  mathematical  statement  of 
the  inverse  problem  and  represent  it  as  a  scattering  problem 
for  the  SchrOdinqer  equation.  We  conclude  the  chapter  with  a 
review  of  the  relevant  literature. 

The  input  to  the  inverse  procedure,  the  plane-wave  re¬ 
flection  coefficient,  and  particularly  its  asymptotic 
behavior  for  large  vertical  wavenumbers  are  the  subject  of 
Chapter  III. 

Chapter  IV  presents  an  extension  of  the  Gelf and-Levitan 
inversion  method  to  the  case  of  non-zero  terminal  potential. 
It  is  this^  formulation  that  permits  us  an  exact  solution  of 
the  inverse  problem  so  that  what  remains  is  the  numerical 
solution  of  the  integral  equation  characterizing  the 
solution. 

The  derivation  presented  in  Chapter  IV  is  followed  in 
Chapter  V  by  a  discussion  of  three  numerical  methods  to  solve 
the  Gelf  and-Levitan  inteqral  equation:  A.n  improved  Born 
approximation  and  two  finite-difference  methods. 

The  numerical  methods  outlined  in  Chapter  V  were  tested 
on  various  postulated  acoustic  profiles  using  synthetically 
generated  reflection  coefficients.  The  representative  numer¬ 
ical  results,  the  impact  of  sampling,  finite  angle  aperture, 
density,  loss  and  noise  are  discussed  in  Chapter  VI. 

Chapter  VII  comprises  the  conclusion  and  suggestions  for 
future  work. 
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CHAPTER  II 


PROBLEM  FORMULATION  AND  REVIEW  OF  PAST  WORK 


The  determination  of  the  acoustic  properties  of  the 
ocean  bottom  from  the  monochromatic  plane  wave  reflection 
coefficient  at  all  anqles  of  incidence  is  now  shown  to  be 
related  to  a  class  of  inverse  problems  in  quantum  scatterinq 
theory  where  an  unknown  potential  in  Schrbdinqer ' s  equation 
is  souqht  from  scatterinq  data.  The  first  part  of  the 
chapter  casts  the  problem  into  mathematical  form  based  on  the 
assumptions  set  forth  in  Chapter  I.  This  is  followed  by  a 
review  of  the  relevant  inverse  problem  literature. 

2 . 1  Problem  Formulation 

( a )  Acoustic  Wave  Equation 

In  acoustics,  the  pressure  qradient  qives  rise  to  an 
acceleration  of  mass  density  p  accordinq  to 

P  §|  =  -  VP  (2.1) 

-► 

where  p  is  acoustic  pressure  and  v  is  particle  velocity. 

Mass  conservation,  toqether  with  a  constitutive  relation 
(Hooke's  law),  yields: 

||  =  -  pc2V • v  (2.2) 

in  which  c  is  sound  velocity  in  the  medium. 
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Fourier  transf orminq  the  time  dependence  in  the  funda¬ 
mental  equations  (2.1)  and  (2.2),  [|^-  -►  -  iu>],  and 

combininq  the  resultinq  time-independent  equations  leads  to 
the  acoustic  wave  equation: 

P(z)V.  (— j-j-y  Vp(  x ,  z )  )  +  kQ2n2(z)  p(x,z)  =  0  (2.3) 

c 

with  k  =  —  and  index  of  refraction  n(z)  =  — r^-r  •  1°  the 

o  c  c  ( z ) 

o 

derivation  of  equation  (2.3),  the  acoustic  medium  has  been 
assumed  to  be  vertically  inhomoqeneous  (or  horizontally 
stratified).  In  other  words,  the  material  properties  are  a 
function  of  depth  (z)  only. 

The  neglect  of  density  variations  reduces  equation  (2.3) 
to  the  Helmholtz  equation  for  the  pressure, 

V2p(x,z)  +  kQ2n2(z)  p(x,z)  =  0  .  (2.4) 

The  equation  (2.4)  constitutes  the  startinq  point  of  this 
study . 

Note :  In  the  presence  of  smooth  density  variations,  the 

acoustic  equation  (2.3)  can  also  be  reduced  to  the 
Helmholtz  equation  through  the  chanqe  of  variable^*^ 


-17- 


P  =  P/Vp 


It  follows  that. 


V2P  +  k  2n,2(z)P  =  0  (2.5) 

o 

with 

,  2  2  -2  ,1  2  3  ,  1  v  2> 

n  =  n  +  ko  C^-p  *  p  “  J  (p  7p)  )  * 

( b )  Mapping  the  Seabed  Below  a  Homogeneous  Ocean 

The  specific  problem  of  interest  is  mapping  the  seabed 
below  a  homogenous  ocean.  The  starting  point  is  again  the 
acoustic  wave  eguation  (2.4): 

V2p(x,z)  +  kQ2n2(z)  p(x,z)  =  0  (2.6) 

for  the  configuration  shown  in  Fig.  2.  Let 

n2(z)  =  1  +  p2(z) 

Since  the  medium  is  homogeneous  in  x,  the  spatial 
variables  can  be  separated  by  assuming  that: 

i(k  sine)x 

p(x,z,k)  =  u(z,k)  e  (2.7) 

Note  that  kQ  sine  =  kx,  the  horizontal  wavenumber  (cf. 

Fig .  3)  . 
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Substitutinq  from  (2.7)  for  o(x,z)  into  wave  equation 
(2.6),  one  obtains  an  equation  for  the  z  variation  of 
pressure ,  u ( z ) : 

r~iL_  +  ( k  2  -  k  2  sin2e)  +  k  2  y2(z)]u  =  0  (2 

l  , 2  v  o  o  ‘  o  J 

dz 


or 

— +  [kQ2  cos28  +  kQ2  y2(z)]u  =  0  (2.9) 

dz 

this  is  similar  to  the  time-independent  Schrodinqer  equation. 
Equation  (2.9)  can  be  written  in  the  familiar  form: 

(  z ,  E )  +  [E  -  V(  z  )  ]  u(z,E)  =  0  (2.10) 

dz2 

with  the  identifications: 

kz  =  kocos0,  the  vertical  wavenumber  is  "momentum". 

E  =  ko2cos20  is  the  "enerqy" . 

V(z)  =  k02jJ2(z)  is  the  "potential". 

A  few  remarks  are  in  order: 

•  The  vertical  wavenumber,  kz,  becomes  imaginary 

for  kx  >  kQ. 
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•  As  k?  ranqes  on  [kQ,  » ]  ,  kx  becomes  pure 
imao i na  rv . 

•  The  enerqy  E  =  k02C0S2Q  j_s  aiWays  real.  For  kz 
real,  E  is  positive. 

•  The  potential  V(z)  =  <0?(1  -  n2(z))  is  in 
qeneral  positive  except,  possibly,  for  a  low 
velocity  layer  at  the  ocean  bottom  interface. 

By  analoqy  with  quantum  mechanics,  one  can  draw  a  "potential 
well"  diaqram  (cf.  Fiq.  4). 

The  SchrOdinqer  equation  (eq.  2.10)  has  associated  with 
it  two  asymptotic  boundary  conditions: 

ik  z  -ik  z 

u(z)  ~  - —  +  r ( p ) — — - —  z  -*•  -® 

/  2ir  /  2tt 


(2.11) 


ik  7, 

u(z)  ~  t(p)  - — z  +  = 

✓77 

where  r(kz)  and  t(kz)  are  identified  as  reflection  and  trans¬ 
mission  coefficients  respectively. 

The  problem  that  was  proposed  in  the  introduction  has 
now  been  cast  into  an  equivalent  quantum  mechanical  problem: 
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Potential 


Fig.4  Scattering  Potential-  Energy  Diagram 

<E],  E2  represent  possible  Energy  levels) 
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Given  r(kz),  the  reflection  coefficient,  as  a  function 


of  vertical  wavenumber  kz,  obtain  the  scattering  potential 
V(z)  in  equation  (2.10). 

2 . 2  Survey  of  Inverse  Methods 
2.2.1  Introduction 

Much  of  the  background  methodology  relevant  to  our 
problem  is  found  in  the  geophysical  literature.  An  excellent 
review  of  the  field  is  provided  by  Newton^2  *  ^  . 

The  seismic  inverse  problem  for  horizontal  layered  media 
of  infinite  depth  consists  in  determining  the  vertical 
structure  of  the  acoustic  medium  (specified  usually  by 
impedance,  or,  in  more  detail,  bv  density  and  velocity)  from 
reflection  measurements.  But  for  a  few  exceptions,  most  of 
the  previous  analyses  have  been  confined  to  excitation  with 
an  impulsive  pressure  siqnal  (6-function)  and  probing  at 
normal  incidence.  The  Fourier  transform  of  a  6 -function  is 
essentially  flat  in  frequency  domain.  What  is  observed  with 
such  an  excitation  is,  therefore,  the  time  trace  of  the 
resulting  medium  response  or  its  Fourier  transform.  Because 
of  the  assumed  horizontally  layered  structure  of  the  medium, 
and  the  vertical  direction  of  the  signal,  acoustic  properties 
of  the  medium  chanqe  only  with  depth  and  thus  the  problem  is 
one -dimensional. 
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In  what  follows,  we  present  a  review  of  that  past 
inverse  problem  work  that  is  relevant  in  several  different 
respects.  First,  the  discussion  of  the  Gelf and-Levitan  ap¬ 
proach  is  particularly  relevant  because  this  is  our  basic 
approach  in  this  thesis  adapted  to  the  Frisk  experiment. 

Next,  the  discussion  of  prior  work  on  sinqle  frequency 
excitation,  launched  at  non-normal  incidence  presents  the 
state  of  the  problem  before  we  addressed  it. 

The  Deift  and  Trubowitz  method  is  described  because 
Stickler,  havinq  been  briefed  on  our  work,  adapted  the  Deift 
and  Trubowitz  method  to  the  Frisk  experiment  and  was  able  to 
devise  an  alternate  approach  to  its  analysis. 

The  discussion  of  the  Schur  algorithm  reviews  the 
analysis  by  Yaqle  and  Levy  of  probing  with  an  impulsive 
excitation  also  at  non-normal  incidence.  We  comment  on  why 
Yaale  and  Levy  dismissed  the  Gelf and-Levitan  approach  as 
inferior  to  the  Schur  algorithm  although  we  have  in  fact, 
successfully  adapted  Gelf and-Levitan  to  the  solution  of  the 
monochromatic,  non-normal  incidence  problem. 

The  Riccati  equation  method  is  discussed  althouqh  it  was 
not  used.  We  considered  this  approach  and  believe  it  to  be 
promising,  but  this  method  was  not  fully  explored. 
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2.2.2  The  Gelfand-Levitan  Approach 

The  problem  that  was  proposed  at  the  end  of  Section 
(2.1),  the  reconstruction  of  the  potential  of  a  Schrbdinqer 
equation  from  the  reflection  coefficient,  is  related  to  a 
classic  problem  of  quantum  scatterinq  theory:  How  to 
reconstruct  a  Sturm-Liouvi 1 le  differential  equation  from  its 
spectral  function.  The  problem  was  solved  in  a  celebrated 
paper  by  Gelfand  and  Levitan  (1^)  w|-10/  since  they  were 
discussinq  the  radial  wave  equation,  were  interested  only  in 
a  solution  on  the  half-line  0  <  r  <  ®  ( stand inq-wave 
problem).  Subsequent  developments  (e.q.,  (15))  led  to 
formulations  on  the  full  line  -«  <  z  <  ®  in  terms  of  such 
readily  measured  quantities  as  the  phase  shift  or  reflection 
coefficient  ( travelinq-wave  problem).  An  excellent  distil¬ 
lation  of  these  ideas  is  to  be  found  in  the  papers  of 
Faddeyev^ ^ ,  while  a  more  qeneral  survey  of  the  field  of 
inverse  scatterinq  has  been  carried  out  more  recently  by 
Chadan  and  Sabatier^7).  The  interrelation  between  the  dif¬ 
ferent  approaches  and  their  time-domain  interpretation  has 
been  presented  by  Burridqe^^.  A  detailed  theoretical 
presentation  and  extension  of  the  Gelfand-Levitan  theory  and 
its  application  to  our  problem  will  be  taken  up  in 
Chapter  IV. 

The  exploitation  of  the  Gelfand-Levitan  formalism  out¬ 
side  of  quantum  mechanics  was  first  taken  up  by  Kay^l^)  and 
Moses  and  deRidder^O)  to  solve  problems  in  electromaqnetics 
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such  as  the  characterization  of  transmission  lines  and 
dielectrics  from  scatterinq  data.  The  possibility  of  mono¬ 
chromatic  nrobinq  is  mentioned  briefly  but  not  pursued.  An 
interestinq  by-product  of  this  research  is  the  theoretical 
construction  of  dielectrics  which  are  reflectionless  at  all 
f  requencies . 

The  Gelf and-Levitan  approach  was  introduced  into  the 
field  of  seismic  exploration  by  Ware  and  Aki^^.  They 
presented  an  analytic  approach  to  the  inverse  scatterinq 
problem  for  elastic  wave  propaqation  in  a  stratified  medium 
when  the  medium  is  probed  with  impulsive  plane  waves  at 
normal  incidence.  The  analytic  solution  was  obtained  by 
transforming  the  equation  of  motion  in  a  stratified  elastic 
medium  for  plane  waves  at  normal  incidence  into  a  one-dimen¬ 
sional  Schrodinqer  equation.  The  potential  of  the  resulting 
Schrodinqer  equation  depends  only  on  the  impedance  of  the 
medium  as  a  function  of  travel  time.  No  ambiguity  arises 
owing  to  the  bound  state  solutions  of  the  Schrodinqer  equa¬ 
tion.  Ware  and  Aki  went  on  to  establish  a  discrete  analoqy 

/ 

of  the  continuous  solution  showinq  aqain  that  the  impedance 
of  the  medium  could  be  recovered  as  a  result  of  probing  at 
normal  incidence  when  the  medium  consists  of  a  homogeneous 
half-space  of  impedance  p0cq  in  contact  with  a  sequence  of 
n  homogeneous  layers  of  impedance  P  jc^  r p  2C2  •  •  •  •  *0  ncn  an<^ 
terminates  with  a  homogeneous  half-space  of  impedance 
pn+lcn+l*  T^e  sec*uence  of  n  homogeneous  layers  which  have 
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thickness  An^  and  which  are  chosen  such  that  travel  time 
throuqh  each  layer  is  a  constant, 

A  n i /c i  =  At  =  constant 

constitutes  what  is  known  as  a  Goupillaud  layered  medium. 

An  eleqant  solution  or  the  inversion  of  a  Goupillaud 
medium  has  been  given  by  Claerbout ^ ^2 )  usinq  z-transf  orms . 
Ware  and  Aki  showed  the  equivalence  of  the  Goupillaud  solu¬ 
tion  and  of  the  discretized  version  of  their  continuous 
solution.  They  had  promised  a  second  paper  dealinq  with  the 
inverse  scatterinq  problem  for  plane  waves  at  non-normal 
incidence.  Such  a  paper  was,  however,  never  published. 
Althouqh  Ware  had  obtained  in  his  thesis  some  partial  results 
at  non-normal  incidence  prior  to  the  publication  of  the  Ware 
and  Aki  paper,  the  aoproach  in  the  thesis  was  too  cumbersome 
and  in  fact  had  hit  an  unsurmountable  wall  at  and  above  the 
critical  anqle:  the  reflection  coefficient  tends  to  one  as 
o,  +  «#  and  therefore  fails  to  meet  an  inteqrability  criterion 
required  in  the  application  of  the  Gelf and-Levitan  algorithm. 

The  question  arises  as  to  how  our  approach,  a  monochro¬ 
matic  experiment  at  all  anqles  of  incidence,  relates  to  the 
Ware  and  Aki  experiment  of  an  impulsive  broadband  source  at 
one  anqle  of  incidence?  The  vertical  wavenumbers  qenerated 

in  Ware  and  Aki  (normal  incidence)  k  =  ~  (for  0  <  u>  <  <*> )  • 

z  co 

Cover  the  ranqe  from  0  to  »  as  the  frequency  is  swept.  One 
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can  see  from  Fiq.  5  that  in  our  experiment  k,  =  k^cose  (for 
0  <  cose  <  »)  formally  covers  the  same  ranqe  of  values  of  kz 
althouqh  the  experiment  is  monochromatic.  There  is,  however, 
a  fundamental  difference  between  the  two  approaches  in  that, 
as  we  prove  in  Chapter  III,  the  reflection  response  in  our 
monochromatic  experiment  is  inteqrable  and  in  fact  qoes  to 
zero  as  kz  qoes  to  infinity  at  least  as  fast  as  (l/kz2).  It 
should  be  noted  that  Ware  and  Aki  did  not  run  any  computer 
simulations  of  their  alqorithm,  and  were  therefore  unaware  of 
its  numerical  performance  (in  fact,  the  Gelf and-Levitai, 
approach  was  widely  held  at  the  time  to  be  numerically 
unstable ) . 

Inspired  by  the  Ware  and  Aki  approach,  a  number  of 
researchers  particularly  Ahn,  Jordan,  and  Krit ikos ^ 2 3 '24 '2 5 ^ 
applied  the  Gelf and-Levitan  alqorithm  to  the  analytical 
problem  of  the  reconstruction  of  dielectric  functions  and 
electron  density  profiles.  Their  work  is  an  analytical 
attempt  to  solve  the  problem  when  a  dielectric  medium  is 
probed  with  impulses  at  normal  incidence.  Most  of  their 
effort  was  applied  to  the  closed-form  solution  of  the 
Gelf and-Levitan  equation.  Such  a  solution  is  possible  when 
the  reflection  coefficient  can  be  represented  as  a  rational 
function  of  wavenumber.  Althouqh  the  approximation  of  the 
reflection  coefficient  by  rational  functions  has  not  yet 
received  any  practical  application,  the  availability  of  such 
closed-form  solutions  provided  us  with  valuable  canonical 
examples  aqainst  which  to  check  numerical  inversion  results. 
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Fig.  5  Generation  of  Vertical  Wavenumbers  in  the 

Ware  and  Aki  Method  and  in  the  Frisk  Experiment 
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More  recently,  Berryman  and  Greene^6)  ^ave  addressed 
some  lonq-standinq  questions  reqardinq  the  qeneral  applica¬ 
bility  of  the  Goupillaud  method.  They  demonstrated  the 
equivalence  of  the  Goupillaud  method  of  inversion  and  of  the 
Marchenko  method^’’)  for  the  SchrOdinqer  equation  for  models 
with  arbitrary  layer  thicknesses  (i.e.,  continuous  impedance 
variation).  When  the  reflection  coefficients  are  correctly 
interpreted,  in  the  continuum  limit,  both  methods  will 
reconstruct  the  same  impedance  except,  possibly,  for  the 
values  at  a  finite  number  of  jump  points  in  any  finite  span 
of  travel  time.  As  part  of  this  work,  Berryman  and  Greene 
presented  a  fast  (0(N2))  recursive  alqorithm  analoqous  to  the 
Levinson  procedure  for  the  inversion  of  a  Toeplitz  matrix. 

We  were  able  to  adapt  this  alqorithm  and  use  it  in  our  numer¬ 
ical  computations. 

Durinq  our  research,  we  became  aware  of  an  unpublished 
report  by  Jacobs  and  Stolt^27^  which  demonstrates  four 
different  coordinate  transformations  which  convert  the 
laterally  homoqeneous  acoustic  wave  equation  of  the 
SchrOdinqer  form.  One  of  the  transformations  takes  frequency 
iii  to  be  a  fixed  parameter  which  infers  our  monochromatic 
condition.  However,  Jacobs  and  Stolt  use  a  sliqhtly 
different  potential  function  than  the  one  chosen  in  this 
thesis.  Their  effort  to  verify  the  Gelf and-Levitan  alqorithm 
is  similar  to  the  one  presented  earlier  by  Moses  and 
deRidder ( 20 )  .  in  discussinq  the  Gelf and-Levitan  alqorithm, 
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they  make  the  important  assertion  that  the  alqorithm  holds 
even  for  the  case  of  dissimilar  end  potentials.  We  have 
presented,  independently,  a  riqorous  proof  of  that  assertion 
in  Chapter  IV. 

The  method  of  Carroll  and  Santosa^8)  was  use(j  by 
Santosa ^ ^ ^  to  solve  the  inverse  problem  for  an  impulsive 
source  at  normal  incidence.  Althouqh  the  method  is  similar 
to  the  Gelfand-Levitan  approach,  it  is  not  based  on  the 
Schrodinqer  equation  but  rather  on  the  equation 


A 


2  A 
U)  v 


A 

q(y ) v ' 


(2.12) 


where  v  is  the  Fourier  transform  of  the  shear  displacement 
and  q(y)  is  related  to  the  chanqe  in  the  impedance  A(y)  by 


q(y)  = 


A1  (y) 
Alyl 


The  measured  response  q(t) 


(2.13) 

6(t,  x  =  0)  is  transformed 


into  the  spectral  density  G(oj) 


The  Gelf and-Levitan  type  equation  that  is  to  be  solved  is 
then 


y 

T(y,z)  +  K(y,z)  -  /  K(y,n)T  (n»z)dn  =  0  ,  z  <  y 

0  n 


(2.15) 


where  T(y,z)  is  given  by 


T(y,z)  =  / 
0 


sin  a)Z 
& j 


cosajy[G((o)  -  — ] du 


(2.16) 


and  the  impedance  profile  is  recovered  from  K(y,x)  through 


a(y) 


2[K(y,y.)l ; 
[1  -  K(y,y) ] 


(2.17) 


The  major  difference  from  the  standard  Gelf and-Levitan 
procedure  is  that  in  (2.15)  the  kernel  T  is  differentiated 
with  respect  to  n . 

Santosa^0^  refined  the  method  to  give  it  a  time-domain 
meaninq  by  applying  it  to  problems  in  which  the  response  data 
are  given  for  a  finite  time.  The  representation  obtained  is 
similar  to  that  in  the  Gopinath-Sondhi  equation^^.  Santosa 
demonstrated  the  method  to  be  stable  both  theoretically  and 
numerically  on  a  window  type  profile.  Reconstruction  errors 
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at  depth  are  attributed  to  errors  in  the  reflection  data  and 
to  the  first  order  discretization  errors  committed  in  the 
approximation  of  an  inteqral  by  a  sum. 

Coen  has  extended  in  one  of  his  oapers^^  the  work  of 
Ware  and  Aki  so  as  to  recover  both  the  density  and  compressi¬ 
bility  profiles  of  a  layered  fluid  from  the  plane  wave 
reflection  coefficient  at  two  precritical  angles  of 
incidence,  and  at  all  frequencies. 

In  another  paper(^3)^  Coen  applied  the  Ware  and  Aki 
method  to  recover  the  three  elastic  profiles  of  a  layered 
half-space  from  three  reflection  coefficients.  First  the 
shear  modulus  and  density  profiles  are  determined  from 
reflection  coefficient  data  for  oblique  incidence  SH  plane 
waves  qiven  at  two  anqles  of  incidence  and  for  all 
frequencies.  Once  the  density  and  shear  modulus  have  been 
obtained,  a  further  experiment  usinq  the  reflection  coeffi¬ 
cient  due  to  an  impulsive  normally  incident  P-wave  permits 
the  retrieval  of  the  P-wave  velocity  and  hence  of  the  Lame 
profile.  The  limitation  in  Coen's  work,  as  in  Ware  and 
Aki's,  is  that  the  potential  V(^)  satisfy  the  inteqrability 
condition . 


/  (1  +  |  C  |  )  |  V(  5  )  |  dj;  <  »  (2.18) 

0 
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which  would  be  very  restrictive  in  practice.  One  result  we 
demonstrate  in  our  work  is  that  the  Oelf and-Levitan  algorithm 
still  applies  to  the  fundamental  case  when  V(^)  tends  to  a 
non-zero  finite  value  at  infinity  in  violation  of  the  above 
inteorability  condition. 

In  a  third  paper^^,  Coen  addresses  the  problem  of 
common  source  point  surface  data  wherein  a  source  is  placed 
on  the  free  surface  of  a  plane  stratified  half-space  and  the 
vertical  component  of  velocity  or  of  acceleration  is  measured 
on  the  free  surface.  After  solvinq  the  impulsive  source 
problem,  Coen  discusses  the  monochromatic  source  problem. 

His  approach  is  deceptively  similar  to  the  one  we  present  in 
this  thesis  as  both  approaches  transform  the  original  problem 
into  a  one-dimensional  SchrOdinger  equation  and  then  proceed 
to  use  the  Gelf and-Levitan  inteqral  equation  to  solve  the 
inverse  problem.  However,  the  problem  in  the  two  approaches 
is  posed  in  a  different  way  and  the  steps  towards  the  solu¬ 
tion  are  dissimilar. 

It  is  useful  here  to  run  through  Coen’s  method  so  as  to 
point  out  the  difficulty  he  encounters  and  which  does  not 
arise  in  our  approach  (a  constant  density  is  assumed). 
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The  Hankel  transform  of  the  pressure  field 


$  (  C  /  z  /  u 1) 


/  f  p(r,z,t)rJn(r/5  )e”*utdrdt 
0  0 


(2.19) 


satisfies  the  Schrodinqer  equation 


jv  ^  O 

( — 5-  -  k  )$(k,z,w)  =  0(  z,u»  )♦  (k  ,  z  ,u  )  .  (2.20) 

3  z 


where  Q(z,uj)  is  the  potential 


0(  z  ,  a) )  = 


^(1 


(2.21) 


and  k  is  related  to  the  horizontal  wavenumber  5  throuqh 


<  =  /k  +  T- 


(2.22) 
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The  SchrOdinqer  equation  is  accompanied  by  two  initial 
conditions 


a  A  =  1 

*  2=0 


z=0 


q  U  #~u5 T 


(2.23) 


where 


q(5'“)  =  ~  prffTsrr.'ST 


(2.24) 


and  d(c,u>)  is  the  Hankel  transform  of  the  vertical  component 
of  particle  acceleration  at  the  surface  z  =  0. 

Coen's  scheme  proceeds  from  an  input  function  r(k,o>) 


r( k ) 


_  1  -  kq(kyU) 
1  +  kq  ( k  ,“(JT 


(2.25) 


qiven  for  all  real  positive  k  values  and  requires  the  compu¬ 
tation  or  R(z,w)  where 


r( k ,u) 


R ( z  ,oj  )e“kzdz 


(2.26) 
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The  Gelf and-Levitan  equation 


z 

A  (  z  ,  y  ,  o) )  =  R(z+y,co)  +  f  A( z , x ,w )R( y+x ) dy  |y|  <  z 

-z 


(2.27) 

is  then  solved  for  A(z,y)  which  in  turn  yields  the  potential 

0(z)  =  2  A(z,z)  ,  z  >  0  .  (2.28) 

The  difficulty  with  the  whole  procedure  stems  from  the 
second  step,  namely  the  computation  of  R(z)  from  r(k).  That 
involves  an  inverse  Laplace  transform  which  is  numerically 
inherently  unstable.  Our  method,  on  the  other  hand,  starts 
off  from  the  Schrodinqer  equation  for  the  field  (rather  than 
for  its  Hankel  transform)  with  the  associated  plane  wave 
reflection  coefficient  as  a  function  of  vertical  wavenumber. 
The  Laplace  transform  of  Coen's  approach  is  replaced  by  a 
Fourier  transform  which  does  not  present  any  numerical  diffi¬ 
culties.  It  is  to  be  noted  that  the  known  numerical 
instability  of  the  Laplace  transform  has  led  some 
researchers,  notably  Santosa  and  Symes(42)  to  dismiss  the 
Gelf and-Levitan  approach  to  the  solution  of  the  inverse 
problem.  We  believe  that  our  approach  to  the  inverse  problem 
could  lead  to  a  positive  reassessment  of  the  Gelf and-Levitan 
inverse  method. 
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2.2.3  Single  Frequency,  Non-Normal  Incidence 

Very  few  researchers  other  than  Coen  have  considered  the 
single  frequency  non-normal  incidence  case.  One  exception  is 
Mittra  and  Schaubert ^ 33 ^  who  used  a  method  different  from 
ours.  Their  approach  is  a  spectral  domain  method  of  probing 
stratified,  lossless,  dielectric  media  using  an  alternative 
to  the  Marchenko  formulation  and  resulting  in  a  Fredholm 
equation  of  the  second  kind  which  is  solved  through  the  use 
of  rational  basis  functions.  They  noted  that  accurate 
inversions  can  be  obtained  if  data  is  provided  for  k,  >>  k  . 
The  Mittra  and  Schaubert  examples  all  have  zero  terminal 
potentials,  and  although  the  results  are  good  in  general,  the 
inaccuracies  are  interestingly  larger  near  the  origin  with 
higher  frequencies  "seeming  to  give  better  resolution. " ( 3^ ) 

Another  example  of  the  prior  single  frequency,  non¬ 
normal  incidence  analysis  is  provided  by  the  work  of 
Roger^36^.  Roger  sought  to  determine  the  index  profile  of  a 
dielectric  plate  backed  by  a  perfectly  conducting  plane. 

That  last  fact  complicates  the  problem,  since  the  potential 
is  always  negative  and  bound  states  due  to  surface  waves 
might  exist.  Roger  starts  from  a  nonlinear  integral  equation 
which  he  linearizes  to  obtain  a  Fredholm  equation  of  the 
first  kind  whose  solution  constitutes  an  ill-posed  problem 
(in  the  sense  of  Hadamard).  Roger  solves  this  equation  by 
using  the  Tikhonov  regularization  method.  The  method  fails 
when  the  permittivity  e (z )  exceeds  a  constant  by  more  than 
20%  and  also  when  the  layer  is  thicker  than  1 . 5x . 
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2.2.4  The  Deift  and  Trubowitz  Method 


Deift  and  Trubowitz^7)  introduced  the  trace  method  for 
determining  the  potential  in  one-dimensional  scattering 
problems  for  the  SchrOdinger  equation.  The  trace  method 
requires  as  data  the  reflection  coefficient,  while  the  method 
we  adopted  requires  the  Fourier  transform  of  the  reflection 
coefficient.  Furthermore,  the  trace  method  requires  the 
solution  of  a  nonlinear  differential  equation  while  the 
Gelf and-Levitan  (or  Marchenko)  equation  that  we  use  is  a 
linear  integral  equation. 

Stickler  visited  us  in  Woods  Hole  and  became  interested 
in  adapting  the  Deift  and  Trubowitz  trace  method  to  our 
problem.  He  took  the  same  input,  i.e.,  the  measurement  of 
the  pressure  field  as  a  function  of  range,  where  both  the 
real  and  imaginary  parts  of  the  pressure  field  are  needed. 
After  the  reflection  coefficient  R(k)  is  derived  by  using  the 
same  approach  as  ours,  he  then  introduces  an  auxiliary 
potential,  q(z),  which  is  determined  by  using  the  trace 
formula  methods  of  Deift  and  Trubowitz.  Stickler^7)  defines 
the  auxiliary  potential,  q(z)  by 

_  • 

q (z  )  =  q(z)  -  ^  (-^)  (2.29) 

P 

where  primes  denote  derivatives  with  respect  to  z.  The  auxi¬ 
liary  potential,  q(z),  can  be  determined  from  the  Deift  and 
Trubowitz  trace  formula 
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(2.30) 


q(z)  =  -  —  /  klm[ R(k )u~^(z ,-k) ] dk 

it  o  * 

The  Jost  function  U2<z,k)  in  (2.30)  is  determined  by  solving 

-u2  +  q(z)  u2  =  k^u2  (2.31) 

with  the  boundary  condition 

u2  ~  e  1*cz ,  z  ♦  .  (2.32) 

Deift  and  Trubowitz  have  shown  that  the  usual  iteration 
scheme  for  solving  two  coupled  nonlinear  integral  equations 
such  as  (2.13)  and  (2.14)  converges.  In  our  case,  instead  of 
(2.13)  and  (2.14),  we  solve  for  q(z)  using  a  Gelf and-Levitan 
linear  integral  equation. 

Stickler  presented  two  numerical  examples^  of  applying 
the  Deift  and  Trubowitz  algorithm  on  a  twice  continuous 
function  (in  Chapter  VI  we  apply  our  method  to  one  of  his 
examples  and  refer  to  it  as  the  "Stickler's  Profile").  Since 
Stickler  generated  the  reflection  coefficient  from  the  solu¬ 
tion  of  a  Riccati  equation,  he  had  control  over  the  local 
tolerance  for  the  determination  of  the  reflection  coeffi¬ 
cient.  As  in  the  Gelfand-Levitan  method,  the  results  are 
excellent  for  z/L  <<  1,  but  deteriorate  gradually  with  depth 
(z/L  >  I).  Stickler  attributes  the  degradation  of  his  method 
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to  the  lack  of  accuracy  of  the  local  reflection  coefficients 
generated  via  the  Riccati  equation. 

When  we  compare  Stickler's  numerical  results  to  ours,  we 
observe  the  general  similarity  of  his  results  to  the  ones 
described  in  this  thesis.  This  similarity  should  not  be  too 
surprising  in  view  of  the  close  relationship  between  the  two 
methods  (see  further  discussion  in  Chapter  IV).  We  have  no 
data  to  assess  the  computational  efficiency  of  Stickler's 
method  versus  ours. 


2.2.5  Schur  Algorithm 

Yagle  and  Levy^88'89)  have  adopted  an  algorithm  which 
reconstructs  the  unknown  acoustic  medium  layer  by  layer 
(layer  stripping  procedure).  The  method  is  analogous  to  the 
downward  continuation  method,  in  that  successive  up  and  down- 
qoinq  waves  are  measured  at  the  surface.  The  first 
reflection  of  the  impulse  yields  information  about  the  medium 
immediately  beneath  the  surface  (at  depth  A).  This  informa¬ 
tion  is  used  to  update  the  waves  at  depth  A  which  then 
becomes  the  new  reference  surface.  The  procedure  is  succes¬ 
sively  repeated  until  the  depth  of  interest  is  reached.  The 
Schur  algorithm  applies  to  the  study  of  the  two  component 
system  of  coupled  differential  equations 


qlx(x,t)  +  qlfc(x,t)  =  -  r(x)  q2(x,t) 
q,  (x,t)  -  q-.(x,t)  =  -  r ( x )  q,(x,t) 


(2.33) 
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where  the  subscripts  x  and  t  denote  derivatives  with  respect 
to  x  and  t,  and  r(x),  the  reflectivity  function,  provides  a 
coupling  between  the  downgoing  wave  q^(x,t)  and  the  upgoing 
wave  q2 ( t)  (unit  velocity). 

Yagle  and  Levy  begin  their  derivation  with  the  set  of 
equations  arising  after  an  initial  impulse  excitation  6(t), 
so  that  qj_(x,t)  and  q2(x,t)  can  be  written  as 

q]_(x,t)  =  6  ( t-u )  +  q  ^  (  x ,  t )  u_1(t-u) 

(2.34) 

q2(x,t)  =  q2(x,t)  u_1(t-u)  , 

in  which  causality  has  been  used  (no  waves  exist  for  t  <  u). 
From  (2.33)  and  (2.34)  Yagle  and  Levy  derive 

r(x)  =  2q2 (u , x)  (2.35) 

The  equations  (2.33)  and  (2.35)  constitute  the  continuous 
parameter  fast  Cholesky  recursion  where  qj(x,t)  and  q2(x,t) 
are  updated  to  yield  r(x)  from  equation  (2.35). 

At  this  point,  the  application  of  the  Schur  method 
entails  taking  the  Fourier  transforms  of  the  system  in 

A 

(2.33).  Denoting  the  transform  of  q  by  q,  we  get: 

A  A  A 

qlx  =  “  r  (x)  q2(x,a>) 

(2.36) 

A  A  A 

q2x  =  -  c(x)  q^(x,a)  +  iuq2(x,uj) 

Yagle  and  Levy  thus  find  a  reflection  coefficient 
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R(x,w) 


q2 ( ) 


(2.37) 


Q1 (x/U  ) 


which  obeys  a  Riccati  equation 


A 


2iuR(  x ,  a) )  +  r  (  x  ) 


(R(x,u>)2 


1) 


(2.38) 


with 


r(x)  =  lim  [2iw  R(x,<u)]  .  (2.39) 

0)  ♦  » 

Equations  (2.36),  (2.37),  and  (2.38)  constitute  the  Schur 
algorithm,  while  (2.38)  and  (2.39)  represent  a  continuous 
parameter  dynamic  deconvolution  algorithm.  It  is  to  be  noted 
that  the  discretized  Schur  algorithm  is  similar  to  the  fast 
recursion  procedure  of  Berryman  and  Greene. 

Yagle  and  Levy  assert  in  the  concluding  section  of  their 
paper  that  their  Schur  algorithm  is  computationally  superior 
to  the  Gelf and-Levitan  algorithm  as  used  by  Coen.  This 
observation  may  be  true  of  the  Gelf and-Levitan  procedure  for 
impulsive  sources  at  non-normal  incidence  as  presented  by 
Ware  and  Aki^21^  and  by  Coen^33^,  but  it  certainly  does  not 
apply  to  our  approach.  Yagle  and  Levy's  main  objection  to 
Gelf and-Levitan  is  that  the  boundedness  of  the  potential 
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00 


(2.40) 


/  (1  +  | T I ) |v(t) jdT  < 

0 

may  not  be  satisfied  at  non-normal  incidence.  However,  we 
show  in  Chapter  IV,  that  the  Gelfand-Levitan  algorithm  can  in 
fact  be  applied  even  for  a  non-zero  end  potential  which  would 
render  the  integral  in  (2.40)  infinite.  Moreover,  the  Frisk 
experiment  studied  here  is  performed  at  a  single  frequency 
and  not  with  an  impulsive  excitation.  We  have  shown  that  the 
monochromatic  reflection  coefficient  as  a  function  of  verti¬ 
cal  wave  number  is  integrable  so  that  the  negative  comment  of 
Yagle  and  Levy  does  not  apply  to  our  work. 

Although  not  stated  in  their  paper,  the  application  of 
the  Schur  algorithm  to  the  inverse  problem  in  a  layered 
acoustic  medium  involves  implicitly  approximations  similar  to 
those  inherent  in  the  Claerbout's  migration  method.  Although 
the  Schur  algorithm  constitutes  an  improvement  over  migration 
in  so  far  as  the  downgoing  wave  strength  is  modified  by  the 
upgoing  wave  strength  it  is  still  an  approximation.  Indeed, 
in  some  of  our  earlier  unpublished  work,  we  succeeded  in 
improving  the  Claerbout  migration  method  precisely  by 
introducing  the  coupling  between  the  reflected  and  the  down¬ 
going  wave.  By  contrast,  our  Gelfand-Levitan  approach  is  not 
an  approximation  and  except  for  numerical  computation 
represents  an  exact  formulation  of  the  problem. 


-45- 


2.2.6  Riccati  Equation  Method 


One  method  which  we  considered  upon  Stickler's  sugges¬ 
tion  and  which  appears  to  have  considerable  potential  is 
based  on  the  Riccati  equation.  Reflectivity  as  a  function  of 
travel  time  t  obeys  the  Riccati  equation'  ' 

~  =  -  2iur  -  y ( t  ) ( 1  -  r2)  (2.41) 

with  boundary  condition  r(oi,»)  =  0  in  which 

(2.42) 


where  the  acoustic  impedance  has  been  defined  by 


Z(x) 


(2.43) 


The  inverse  problem  is  here  that  of  reconstructing  Z(z) 
from  surface  observations  of  r.  To  get  Z(z)  from  Z(t) 
involves  further  assumptions.  We  can  formulate  the  problem 
equivalently  by  writing 


rU,t)  =  /  dt'  y(t'  )el2io,(T,"T)]  (l-r2(u),T))  (2.44) 

T 

while  at  the  surface,  the  reflection  amplitude  is 


rU,o)  =  /  dTY(T)e(i“T)(l-r2U,T)) 
0 


(2.45) 
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The  last  two  equations  can  be  regarded  as  a  nonlinear 
mapping  of  r(w,0),  given  as  a  function  of  u ,  into  y  (x  ) .  The 
aim  of  such  an  approach  would  be  to  reconstruct  y  (t  )  from 
r(u>,0).  Although  the  existence  of  a  complete  solution  to  the 
above  problem  has  not  yet  been  demonstrated,  an  approximation 
can  be  derived.  An  inverse  Fourier  transform  of  (2.28) 
yields 


1  " 

r(t)  =  —  /  due 

■n  * 


-(2io,T) 


r (u ,0 ) 


+  —  /  du  /  dt  'o  (t  '  )r^(u  ,t  ' )  xe^*w^T  T^. 

IT  r\ 

—08  (J 


(2.46) 


The  first  approximation  to  the  solution  is  the  Born  series 
term,  namely 


Yi(t) 


i  00 
-  /  du 

n  1 


r (u , 0  )e 


~2  iut 


(2.47) 


which  resembles  the  first  term  of  our  own  result.  The 
convergence  of  the  iteration  series  solution  has  not  yet  been 
established,  but  the  approach  as  formulated  by  Nilsen  and 
Gjevik(4°)  appears  promising. 
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It  is  useful  here  to  note  that  a  Riccati  equation  still 


applies  in  the  presence  of  density  discontinuities.  The 
equation  obeyed  by  the  pressure  p  is 

pto'V)'  +  O'2  -  V)p  -  0  ,  <2. 

where  V  is  the  potential. 

The  reflection  coefficient  may  be  obtained  at  the 
surface  from  the  continuity  of  pressure  and  vertical  velo¬ 
city. 


R(k) 


z=0 


p-1p'  -  ikp 
* — ^ 


(2.49) 


The  ratio  of  vertical  velocity  to  pressure,  i.e. ,  the  admit¬ 
tance. 


u 


(2.50) 


is  continuous  even  in  the  presence  of  material  discontinu¬ 
ities,  as  it  is  the  ratio  of  two  continuous  quantities.  The 
admittance  u  obeys  a  Riccati  equation. 
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(2.51) 


in  which  q 


V) 


Note  that  u(z=0)  =  ik-|— 


Let  u  ( z  )  =  i k 

It  follows  that  W(0)  =  R(k)|  0  ,  and  W(z)  itself  obeys  a 

Riccati  equation 


2ikW'  -  2W(0  +  pk2)  +  W'2(pk2  -  0)  =  (0  -  pk2) 

(2.52) 

which  can  be  used  to  generate  the  reflection  coefficient  in 
the  direct  problem  or  to  use  an  iterative  procedure  for  the 
inverse  problem.  We  did  not  pursue  this  approach  further  but 
it  merits  further  study. 

We  have  covered  here  those  papers  that  were  most 
relevant  to  the  work  that  follows.  It  should  be  noted,  how- 
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ever,  that  there  are  a  number  of  other  interesting  inverse 
methods  that  have  not  been  included  in  this  review. 
Particularly  noteworthy  are  the  papers  by  Moses^4^,  Gopinath 
and  Sondhi(3l),  Santosa  and  Symes(42)  and  others ( 4 3 • 4 4 • . 

Also  of  interest  is  additional  literature  on  the  approximate 
inverse  methods ( 55-57 ) .  special  mention  should  be  made  of 
methods  based  on  the  Born  approximation  which  originated  with 
the  seminal  work  of  Cohen  and  Ble istein * 45-54 ) . 
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CHAPTER  III 


THE  REFLECTION  COEFFICIENT 

3 . 1  Introduction 

The  plane-wave  reflection  coefficient  as  a  function  of 
vertical  wavenumber  r(kz)  is  central  to  the  inversion  proce¬ 
dure  we  adopted  to  solve  the  inverse  problem.  Its  symmetry 
property  is  shown  to  follow  from  the  integral  equation  repre¬ 
sentation  of  the  field.  We  then  derive  its  asymptotic  behavior 

for  large  k  by  induction  and  show  that  r(k  )  is  integrable, 
z  z 

and  hence  is  an  acceptable  input  to  the  inverse  method  presented 
in  Chapter  IV. 

3 . 2  Definition  and  Properties  of  the  Reflection  Coefficient 
Consider  the  problem  of  a  plane  wave 


\pi  (z)  =  e 


ik  z 
z 


incident  from  z  =  -»  onto  a  half-space  extending  from  z  =  0 
to  z  =  +°°.  The  half  space  is  characterized  by  a  potential 
V(z).  V(z)  is  all  that  is  needed  to  describe  the  scattering 
of  the  incident  wave  iJk  (z)  by  the  acoustic  half  space: 


ik  z 

(z)  =  e  z  + 


G(z,z')  V(z')  ^(z')dz'  -00  <z  <  00 

(3.1) 
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where  the  Green's  function  G(z,z')  is  determined  by 


G (z , z  ' )  = 


ik.,  z-z' 
e  ^ 


(3.2) 


It  follows  that  as  z 


ik  z  -ik  z  c«>  ik  z' 
iMz)  =  e  2  +  e  2  - — —  V ( z  ' )  ^(z')dz' 

2ik 

2  (3.3) 


ik  z 

and  the  coefficient  of  e  2  can  be  identified  as  the. 
reflection  coefficient. 


1  ( 

r(kz)  =  2ik 

2  -00 


°°  +ik  z ' 

e  2  V (z  1 )  ijj(z  1 ) dz  1 


(3.4) 


The  integral  representation  of  the  reflection  coefficient 

allows  a  simple  derivation  of  the  symmetry  properties  of 

r(k  )  in  the  complex  plane, 
z 

From  (3.4) 


r(-kz)  =  -  jt—  f  e  2  V(z')  ipiz’,  -  dz ' 

2  loo 


(3.5) 


But,  for  a  real  potential  V(z)  and  for  real  k  ,  Schrodinger ' s 

£» 

equation  shows  that  when  ^(z',k  )  is  a  solution,  so  is 

ip  ( z',  -k  ).  Moreover,  ip(  z',  -  k  )  =  ^*(z',k  )  which 

z  z  z 

implies 


52 


(3.6) 


r(-kz)  =  r*(kz) 

The  result  can  be  extended  to  the  complex  plane  through 
the  use  of  the  Schwartz  reflection  principle, 

r(-k  *)  =  r*  (k  )  (3.7) 

in  any  region  of  analyticity  connected  with  the  real  k  axis. 

z 

(This  constitutes  the  analytic  continuation  of  r(k  )). 

z 

Now,  the  Gelf and-Levitan  algorithm  requires 

knowledge  of  r(kz)  for  all  real  kz.  However,  the  symmetry 

property  expressed  in  (equation  3.6)  demonstrates  that 

knowledge  of  r(k  )  on  the  half-line  of  0  <  k  <°°  is 

z  z 

sufficient. 

Note  that  when  the  vertical  wavenumber  k  =  k^cosG 

z  0 

is  real  and  larger  than  the  water  wavenumber  kQ,  the  angle 
of  incidence  becomes  imaginary.  That  can  be  verified  by 
requiring  that  cosO  =  cos(Qr  +  iO^)  be  real.  And  since 

cos(0  +  i0.)  =  cos0  coshG.  -  i  sin0  sinh0. 

i-  1  1>  1  X 

it  follows  that  0^  =  0  and  k  =  kAcosh0.  for  kA  <  k  <  °°. 

r  z  0  i  0  z 

The  mapping  between  the  k  -plane  and  the  0-plane  is  drawn 

Z 

in  Figure  6. 


53 


r r 


Re(kz) 


\r  /  r  r  r  / - r 


Re(8) 


Fig.  6  Mapping  from  the  kz~  Plane  into  the  0-Plane 


3.3 


Asymptotic  Behavior  of  the  Plane  Wave  Reflection 
Coefficient 

The  input  to  both  the  Gelf and-Levitan  algorithm 
and  to  the  Born  approximation  is  the  Fourier  transform  of 
the  reflection  coefficient 


R(z) 


f°°  ik  z 

r(kz)  e  2  dkz 

—CO 


The  properties  of  r  (k  )  as  k  ■+»  are  studied  for  two 

z  z 

simple  cases  and  then  generalized. 


(a)  Half  space 

The  reflection  coefficient  is  given  by 


r(kz) 


where  the  impedance  Z  =  pc/cos 0  or 


(3.8) 


(3.9) 


Therefore, 


r(kz) 


kz  I  kzl 
kz  +  kzl 


(3.10) 
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And  since  is  continuous, 


kzl  =  <kl2  -  +  kz2 


(3.11) 


k  -  V(k  2  -  k  2)  +  k  2' 

r(k  )  =  -2 -  1  °  * 

kz  +  V(kl2  "  k02)  +  kz2 


(3.12) 


As  k  00 
z 


r(kz)  <«, 


kz  -  V1  *  3*7  »ki2  -  ^o2>> 

kz  +  kz(1  +  “^2  (kl2  “  k02) 

2  k 

z 


(3.13) 


=  _  i (fc  2  -  k  2) /k  2 
4  'K1  K0  ,/Kz 


Hence, 


r<V  nr  t  £ 0  as  <!/kz  > 

z 


When  density  variations  are  considered. 


r<kz)  = 


mcosG  -  ncos  0, 


mcosG  +  ncosQ, 


(3.15) 
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where  m  =  p  /p  and  n  =  k^kg 


r(kz) 


-  k 


+  k 


zl 


zl 


As  k  ■+  ® 
z 

r<V  '  +  <kl2  -  k02,/2kz2  (3-17» 

2 

which  still  decays  as  (1/k  ) ,  but  tends  to  a  finite 

z 

limit  <~sir~)  as  kz  *  ” 

Therefore,  the  Fourier  transform  of  r(k  )  involves 

z 

generalized  functions. 


F(r(kz))  =  (j— £y)6(z)  +  (analytic  function)  (3.18) 


(b) 

is 


One  Layer  Case  (cf .  Figure  7) . 

The  reflection  coefficient  at  the  (0  -  1)  interface 


,Z1  I  \  ,  ,Z2  ~  Z1  2iklzd 

'•7  X  9  )  ^  '  T7  X  9  )  ® 


r(kz) 


‘Z1  +  Z0 


Z2  +  Z1 


Z.  -  Zn  z9  -  Z.  2ik.  d 

*  X  17  )  '9  JL  n  )  ^ 


Z1  +  zo 


Z2  +  Z1 


(3.19) 


Fig.  7  One  Layer  Case 
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The  analysis  of  the  terms  within  parentheses  is 
identical  to  the  analysis  carried  out  in  the  previous  section. 
One  can  conclude  immediately  that 


i.  i.  2  xk  d  h 

r(kz)  - >  Nk^  -  kQz)  +  (k  2  -  k.z)  e  2  ] /4k  ^ 

k  -►“> 
z 


0  as  (1/k  ) 

z 


(3.20) 


Note: 

The  delta  function  potential 

V(z )  =  A<$  (z)  (3.21) 


can  be  considered  as  a  limiting  case  of  the  one-layer 
problem.  It  has  associated  with  it  a  reflection  coefficient 


r(kz) 


A 

1  2k  +  iA 
z 


(3.22) 


which  goes  to  zero  as  k  +  ®,  but  only  as  (1/k  )  in  apparent 

z  z 

violation  of  the  result  just  derived.  The  reason  is  that 

in  a  delta  function  potential,  V(z)  ■+<»  which  implies 
2  2 

(k0  -  k^  )  -*•«>,  invalidating  the  binomial  approximation  to 

the  square  root  used  in  section  (a) .  However,  the  formula- 

2 

tion  of  the  acoustic  problem  requires  V(z)  £  kg  which  is 
finite  (  as  long  as  the  frequency  w Q  is  finite) .  Therefore, 
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delta-like  potentials  do  not  arise  and  the  binomial 
approximation  is  valid. 


The  continued  fraction  representation  of  (equation  3.25) 
indicates  that  the  partial  reflection  coefficient 
in  (equation  5.23)  could  represent  the  reflection 
coefficient  due  to  a  complicated  medium  rather  than  to  a 
simple  homogeneous  half  space. 

The  asymptotic  behavior  of  r(kz)  is  deduced  in  a  two 

step  process.  Assume  the  configuration  of  Figures 

with  R.  (k  )  -*■  (— =— )  as  k  +  “ 
l  z  ,2  z 

z 

Adding  a  new  interface  to  the  set-up  (Figure  9  ) 
and  using  equation  (5.24), 

2  42  “i2kzdl 

k  (1  -  k’/B*1)  e  z  1 

k  (k  )  *  +  - 5—5 - 

z  p  k  Z  -2ik  d 

Kz  ^  Z1Kzi  ,  a 

“T"  k  2 

Kz 


:z  8a  2ikzdl 

-f-)  (1  ~  —4  e  z  1 

p  v  H 


a2  2  4ik  d. 

+  6a  e  z  1  + 


(3.26) 


Therefore, 


r(kz)  TT^  0  as  (1/kz  ] 
z 


(3.27) 
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Fig.  8  Stack  of  Isovelocity  Layers 


62 


63 


The  argument  just  presented  is,  in  fact,  the  last  step  in 

a  proof  by  induction  of  the  proposition  that  r(k  )  -*-  0 

z 

2 

as  (1/k  )  in  the  presence  of  a  homogeneous  half  space. 

The  proposition  was  proved  true  for  one  or  two  interfaces 
was  assumed  true  for  an  arbitrary  number  of  interfaces  and 
was  shown  to  hold  for  one  more  interface.  This  result 
should  be  contrasted  with  the  corresponding  situation  in 
the  Ware  and  Aki  experiment;  For  angles  of  incidence 
greater  than  critical,  they  were  confronted  with  the 
fact  that  r(o))  ->1,  as  o'  ■+  <*,  and  could  not  proceed  with  the 
Gelfand-Levitan  inversion  procedure. 
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CHAPTER  IV 


THE  ONE-DIMENSIONAL  INVERSE  PROBLEM 


This  chapter  presents  the  derivation  of  our  approach  to 
the  one-dimensional  inverse  problem.  The  results  yield  the 
methodology  underlying  the  numerical  computations  described 
in  Chapter  V.  The  derivation  exploits  the  equivalence 
between  the  acoustic  problem  and  the  corresponding  quantum 
scattering  problem  as  presented  in  Chapter  II.  A  major  dif¬ 
ference  between  the  two  is  the  boundary  conditions.  Whereas 
in  quantum  mechanics,  the  unknown  slab  is  surrounded  by 
isovelocity  space  (zero  end  potential),  in  the  acoustics 
problem,  as  applied  to  the  ocean  bottom,  differing  velocities 
have  to  be  accommodated  above  and  below  the  slab  (non-zero 
end  potential).  The  solution  to  the  inverse  problem  detailed 
in  this  chapter  consists  of  an  extension  of  Faddeev's 
method^*’)  for  deriving  the  scattering  matrix  in  the  case  of 
zero  end  potential  (or  integrable  potential)  to  the  case  of 
non-zero  end  potential  which  is  representative  of  the  reflec¬ 
tions  from  the  bottom  of  the  ocean.  We  are  able  to  solve  the 
problem  analytically,  in  part,  because  of  its  one-dimensional 
modeling.  Our  technique  is  related  to  that  applied  by 
Stickler  to  a  similar  scattering  problem. 
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4.1  Introduction 


The  determination  of  the  potential  V(z)  in  the 
Schrbdinqer  equation 


— %  +  [p2  -  V(z)  ]  <|>  (  z  )  =  0 
dz 


(4.1) 


from  scattering  data  such  as  the  reflection  coefficient 
constitutes  the  one-dimensional  inverse  problem. 

Two  scatterinq  solutions  of  the  Schrftdinqer  equation  are 
defined  by  their  asymptotic  behavior  (p  is  the  vertical  wave- 
number  also  referred  to  as  kz): 


<|» z,p) 


r 


elpz  +  s12(p)  e~lpz  as 


si;L(p)e 


1PZ 


V 


Z  +  -  ® 

z  +  +  • 


(4.2) 


<|»2<  z  *P) 


f  e"ipz  + 


s21(p)elpz  as 


s 2 2 ( p )  e"iPZ 


z  +  » 

z  +  -  » 


(4.3) 


The  matrix  of  coefficients 


11<P> 

s,2(0) 

21,n) 

s22(n) 
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is  known  as  the  S  matrix  of  the  Schrbdinqer  equation.  The 
element  si2(p)  corresponds  to  the  plane  wave  reflection  coef¬ 
ficient  determined  in  the  WHOt  experiment. 

Faddeev^16B^  has  shown  that  s12(p)  may  determine  all  the 
elements  of  the  S  matrix  and  hence  is  sufficient  to  obtain 
the  scatterinq  potential  V<z).  However,  these  results  have 
an  important  practical  restriction,  namely,  that  the  end 
potential  tends  to  zero;  i.e.,  V(z)  0  as  z  -*•  •  They  are 

reviewed  in  the  first  part  of  this  chapter. 

In  the  second  part  of  the  chapter,  we  present  an  exten¬ 
sion  of  the  theory  to  include  the  qeophysically  siqnificant 
case  of  a  finite  end  potential,  V(z)  +  as  z  +  <*>.  Finally, 

an  appropriate  choice  of  source  frequency  is  shown  to  elimi¬ 
nate  trapped  modes. 

4 . 2  Properties  of  the  Solution  of  the  Schrbdinqer 

Equation  (V^O) 

Two  fundamental  solutions  of  the  Schr’ddinqer  equation 
are  introduced. 


/  \  1PZ 

u1(z,p)  ~  e 

as 

z  -*■  +  00 

(4.4) 

,  .  -ipz 

u  2  ( z ,  p )  ~  e 

as 

z  -*■  -  00 

(4.5) 
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The  method  of  variation  of  parameters  leads  to  a 
representation  of  Uj(z,k)  and  u2(z,k)  as  solutions  of 
Volterra  eauations  of  the  second  kind^^ 


u, (z,p) 


=  e 


1PZ 


-  I 


sini? LllZLl  V(z)u1(z*  ,p)dz* 


(4.6) 


u  2  (  z  f  p )  = 


(4.7) 


But,  since  the  SchrOdinqer  equation  is  symmetrical  in  p, 
u1(z,-p)  and  u2(z,-p)  are  also  solutions  of  (eq.  4.1).  The 
solution  pairs  (u1(z,p),  u^Zr-p)]  and  [u2(z,p),  u2(z,-p)] 
are  linearly  independent  since  their  Wronskians  obtained  from 
their  asymptotic  form  (eqs.  4.4,  4.5) 


w(u1,u1*)  =  2ip 

W(u2,u2*)  =  -2ip 


(4.8) 


are  non-zero  for  p  *  0. 

Now,  any  solution  of  the  SchrOdinger  equation  can  be 
written  as  a  linear  combination  of  two  independent  solutions. 
In  particular, 
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(4.9) 


UjUfP)  =  u2(z,p)c22(p)  +  u2(z,-p)c21(p) 

u  2 ( z  #  p )  *  u  ^(z ,p)c  +  u1(z,-p)c12(p)  .  (4.10) 

The  next  section  examines  some  properties  of  the  coefficients 
c^(k)  which  will  later  be  shown  to  be  closely  related  to  the 
elements  of  the  S  matrix. 

4.2.1  Properties  of  the  Coefficients  c^jCp) 

The  coefficients  c^j(p)  can  be  expressed  as  Wronskians 
by  "takinq  Wronskians"  of  both  sides  in  equations  (4.9  and 
4.10).  for  instance,  from  (equation  4.9) 


w(u1(z,p),u2(z,p)}=  w[u1,u1]c11(p) 

+  w[u^(z,p),u^(z,-p)]c^2(p). 


(4.11) 


We  know,  however,  that  in  (4.11) 


w[uj,u^]  =  0  (linearly  dependent  functions) 
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and 


w[u  j  (  z  ,p)  ,u  ^  (z  ,-p)  ]  =  2ip  (cf.,  eq.  4.8)  .  (4.12) 

It  follows  that 

c 12 ( P )  =  2^7  wtui (Z»P) •  u2(Z,p)]  .  (4.13) 

One  can  show  similarly  that 

c21(p)  *  C 12(p)  (4.14) 

cn(p)  ■  -c22(-p)  —  ’5Tlc  W[  u2  (* ,— p)  r  u  x  ( z  ,p )  ]  (4.15) 

Usinq  the  values  of  c^j  expressed  in  this  form,  and  substi- 
tutinq  eq.  4.9  into  eq.  4.10,  one  qets  the  compatibility 
relations , 

cll(p)c22(p)  +  c21(p)c12 =  1 

(4.16) 
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and , 


cl2(p)c22(o)  +  c  (-p)c21(p>  =  0 


(4.17) 


The  asymptotic  form  of  U2(z»k)  in  (ea.  4.7)  for  z  «° , 


u0(z,p)=  e~ipz[l-  J  eip2'v(z'  )Uj(z'  )  dz '  ] 


j^p  2  oo  t 

+  /  e“ipz'v(z' )u2(z' )dz 


(4.18) 


leads  by  comparison  with  (eq.  4.9)  to  the  identifications: 


c12(p)  = 


1  - 


1 

JTp 


,ipz 


V(z 


>u2(z' 


)dz 


(4.19) 


clx (p) 


e"ipz'v(z’ )u2(z' )dz' 


(4.20) 
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Usina  the  previous  estimates  for  U2  in  eq.  4.5  and  the 

Rienann-Lebesque  theorem,  it  can  be  seen,  that  for  larqe  p, 

.  1.  * 

c12(p)-1  is  o(l/p)  and  Cj^(p)  is  o[  — J .  Moreover,  c22(p)  is 
analytic  in  the  upper  half  plane,  does  not  vanish  on  the  real 
p  axis,  and  has  only  a  finite  number  of  simple  zeros  on  the 
imaqinary  axis. 

4.2.2  Properties  of  the  s^j(k)  Coefficients 

One  is  now  ready  to  return  to  the  oriqinal  scatterinq 
problem  and  its  associated  S  matrix. 

ij>]_(z,o),  the  solution  of  the  scattering  problem  of 
interest,  can  be  written  in  terms  of  the  linearly  independent 
solutions  u2(z,p]  and  u2(z,-p) 


^(2,0)  =  u 2 ( z  , ”P )  +  s  12  (p)u2(  z,p) 


(4.21) 


=  slx (p)u j(z ,p) 


(4.22) 


The  order  symbol  o(  )  is  defined  as  follows: 

f  ( e  )  =  o  (q  ( e  )  ]  as  e  ♦  0 

if  lim  =  0 

e  ♦  0  ^T«T 
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Wri.ti.nq  out  u1(zrp)  in  terms  of  u2(z,±p)  as  qiven  in  eq. 
4.10,  one  obtains  the  compatibility  relations 


Sll(p)  =  c21(p) 

C22(p) 

sl2(p)  =  c“TpT 


(4.23) 


Carryinq  out  the  same  operation  on  <J>2(z,k),  one  obtains 


s22!d)  "  c12(D) 

(4.24) 

,  .  C11<D’ 
s2l<p)  -  STTTpT  • 


From  the  asymptotic  behavior  of  the  Jost  functions  c^j 

* 

one  can  deduce  the  asymptotics  of  the  S  matrix 
sll(p)  =  s22<p)  =  1  +  0( 1/d) 

(4.25) 


s 1 2 ( P )  and  s2 i ^ p ^  =  0(l/p) 


The  order  symbol  0(  )  is  defined  as  follows: 

f  (  e  )  =  0(q  (  e  )  ]  as  e  •*■  0 
if  lim  •— 7— -j-  =  A  ,  0  <  I A I  <  <*> 

.  .  0  q(t)  1  1 


73 


It  can  be  seen  from  eq.  4.24, 


4.16  and  eq.  4.17  that  the 


S  matrix  is  unitary. 


(4.26) 


implyinq  conservation  of  enerqy,  and  that  since  s^j(-p)  = 
s*ii(D), 


slls 


21 


S  12s 


22 


=  0 


(4.27) 


The  coefficients  s-^(p)  are  continuous  for  real 
p,  Sj^tp)  beinq  analytic  in  the  upper  half  plane  except  for 
poles  on  the  imaqinary  axis  ( correspondinq  to  the  zeros  of 
cl2(p)).  Conditions  (4.26)  and  (4.27)  allow  one  to  recon¬ 
struct  the  scatterinq  matrix  from  a  knowledge  of  the  reflec¬ 
tion  coefficient  s12(p).  In  what  follows,  sj2(p)  *s 
identified  with  the  plane  wave  reflection  coefficient  r(kz) 
and  s^j^(p)  with  the  transmission  coefficient  t(k2).  Substi¬ 
tuting  c2j  by  (1/t)  and  c22  by  (r/t)  in  (eq.  4.10),  it 
follows  that 


t(p)u1(z,p)  =  r(p)u2(z,p)  +  u  2 ( z , — p )  .  (4.28) 

The  above  equation  is  the  basis  for  the  derivation  of 
the  Gelf and-Levitan  inversion  method  that  has  been  obtained 
for  a  potential  V(z)  ♦  0  as  z  +  ±m»  Such  a  potential  arises 
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in  the  study  of  dielectric  slabs  by  electromagnetic  probing. 
In  the  earth,  however,  the  velocity  c(z)  tends  to  a  value  c^ 
larger  than  the  surface  velocity  Cq.  Correspondingly,  the 
potential  V(z)  tends  to  a  positive  constant  V^, 


Therefore  the  Gelf and-Levitan  inversion  does  not  really 
apply,  and  a  scattering  solution  is  needed  that  allows  a  non¬ 
zero  end  potential.  In  the  next  section,  we  present  our 
approach  to  this  problem. 

4 . 3  Non  Zero  Final  Potential 

To  accommodate  to  a  non-zero  V^,  we  now  consider  the 
fundamental  solution  of  the  Schrbdinger  eguation  ui(ZfP) 
defined  asymptotically  as 

u^Zjp)  -  e*p  z  z  -*•  <»  (4.30) 


where 

P*  =  ~  vi  • 
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We  assume  Imp’>  0,  in  order  to  satisfy  the  radiation  condi¬ 
tion  at  • .  The  Volterra  inteqral  equation  representation  of 
u^(z,p)  in  eq.  4.6  is  modified  to 


Uj(z,o)=e*p  Z-J  r^-- —  -[  V(  z  )-V1]  Uj  (  z  1  ,p)dz  *  (4.31) 


while  the  inteqral  equation  representation  of  u2(z,p)  is 
unchanqed 


u_(z,p)=e  ipz+  f  s*  n-  -z-2 — V(  z  )  u, ( z  '  ,p)dz  ' 


(4.32) 


4.3.1  Asymptotic  Behavior 


The  key  observations  to  be  made  relative  to  eq.  4.31  is 
that  since  uj(z,p)  is  the  solution  of  a  Volterra  equation  of 
the  second  kind  with  square  inteqrable  kernel,  the  method  of 
successive  approximation  will  converqe.  That  observation  has 
in  fact  been  applied  by  D.  Stickler^7  in  connection  with  the 
Deif t-Trubowitz  inversion  procedure?  thus,  we  have 


T  -i  ivi  i 

V(  z  )  -V .  ipz-  ■*— — z-  -yi— 7(V(z' )-Vn  )dz' 

u,(z,p)=  1 - ie  zp  Zlpz  1  +H.O.T.* 

( 2 ip) z  ( Imp  >  0) 


(4.33) 


H.O.T.  stands  for  higher  order  terms 
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from  which  the  transmission  coefficient  t(p)  can  be  deduced 


t(p)=( 1- 


,V1Z0 


V, 


(  2ib) 


') e 


2p 


+2lO 


f  V(z' )dz'+ 


1 

2  ip 


7  (v(z’ )-v1)dz' 

Z0  +H.O.T. 


(4.34) 


where  Zq  is  arbitrary. 
Similarly,  (Imp  >  0) 


V(z) 


u0 (z,p)=( 1- 

(2ip) 


-lpz- 


2lp 


+  H.O. 


(4.35) 


4.3.2  Inversion  Procedure 

We  are  now  ready  to  present  the  inversion  procedure  for 
>  0.  Note  that  the  basic  relation  (eq.  4.28)  still  holds 
for  Vj  >  0, 

t(p)u1(z,p)  =  r(p)u2(z,p)  +  u  2 ( z , ~p )  (4.36) 

Followinq  Faddeev's^^  case  of  V  =  0,  for  the  case  of 
Vx  >  0,  a  function  h(z,p)  =  U2(z,p)e1Pz  is  introduced.  The 
expression  for  U2(Z/b)  (eq.  4.35)  shows  that  h(z,p)-l  is 
analytic  in  the  upper  half  plane  and  Im(p)  >  0  and  ♦  0  as 
|p|  ♦  <*>.  We  thus  obtain  (e  -*•  0  +  ) 
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.  ,  „  _ .  ,  _  1  r  h(z,p* )-l 

h(z'p)_1  -  TUT  J  ' 


(4.37) 


or  usina  (ea.  4.36) 


h  (  z  ,p)  -  1  =  2^-j-  f  do 


.  ®  t(p’ )u. (z,p* )e-ip  z- 

r,w _ ± _ 


p  1  +p+ie 


1  «  r  { p  ’ )u2(z,p* )e  ip  Z 

2 it  i  f  ^  1  ~  ^p 


p '  +p+ie 


(4.38) 


The  first  inteqral  in  (ea.  4.38)  is  zero  since 


t(p)u1 (z,p)e“ipz-l=  - 


V(z)  „ 

- 7  e 

(  2  ip ) z 


Tip  /  V( z ' ) dz ' 


+  H.O.T.  (4.39) 


thus 


,  «  r (o' )u?( z ,p' )e”ip  z 

h(z'D)  -  1  1  '  57T  - dp‘  ’  (4'38a) 

Comparing  egs.  (4.37)  and  (4.38a),  we  can  express  u2(z,p’)  in 
the  Levin  representation^7^, 

UjtZfP)  =  e  *pz  +  J  K(z,z’)e^pz  dz '  (4.40) 

00 

in  which  the  kernel  of  the  inteqral  does  not  depend  on  p; 
i  .  e . , 
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(4.41) 


K(z,y)  =  “  -J7  /  r(p)u2(z , p)e-ipYdp 

—  00 

Insertinq  (eq.  4.40)  into  (eq.  4.41)  results  in  the  Gelfand- 
Levitan  tyoe  equation 

z 

K(z,y)  +  R(z+y)  +  /  R( z 1 +y ) K ( z , z ' ) dz '  =  0  (4.42) 

y  <  z 


where  R(z)  is  the  Fourier  transform  of  r(p) 


R(z)  =  -57  /  r(p)e-10zdp  . 


(4.43) 


When  (eq.  4.41)  is  substituted  into  the  Schrbdinger 
equation,  it  is  found  that  K(z,y)  satisfies  a  partial  differ 
ential  equation 


3  2  j,  a  2  „2 

— j  ~  V(Z)K 

3z  3y 


0 


subject  to  the  boundary  conditions 


(4.44) 


K ( z ,  -«)  =  0  (4.45) 

V(*)  =  2§|  (Z,Z)  . 
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The  solution  of  R(z)  in  eq .  4.43  consists,  in  qeneral, 
of  both  a  continuous  and  a  discrete  part.  The  discrete  part 
applies  in  the  presence  of  trapped  modes.  In  that  case, 
additional  information  is  required  to  construct  R(z) 

l  00  _  •  -ip.  z 

R  ( z )  =  jjj-  /  r(p)e-ipzdp  +  l  rr^e  .  (4.47) 

i 

where  the  p^'s  are  the  poles  of  r(p)  on  the  positive 
imaqinarv  axis. 

The  choice  of  the  constant  normalization  coefficients 
is  dictated  by  the  requirement  that  V(z)  =  0  for  z  <  0.  This 
can  be  seen  by  examininq  the  Gelf and-Levitan  equation  (4.42), 

z 

K(z,y)  +  R(z  +  y)  +  /  R(z'  +  y)  K(z,z')dz'  =  0.  (4.48) 

— oo 

y  <  z 

We  note  that  R(z)  =  0  for  z  <  0  insures  that  K(z,y),  and 
hence  that  V(z)  =  2  (z,z)  are  all  zero  for  z  <  0.  The 

choice  of  m^  in  (eq.  4.47)  is  therefore  dictated  by  the 
requirement  R(z)  =  0  for  z  <  0.  Now,  the  inteqral  in  (eq. 
4.47)  is  for  z  <  0, 

/  r(p)e  *pzdp  =  i  I  b^e  *  (4.49) 

_u>  i 

where  the  b^'s  are  the  residues  at  the  poles  p^  of  r(p). 
Substitutinq  for  the  value  of  the  inteqral  in  (eq.  4.47),  and 
imposinq  r(z)  =  0  for  z  <  0  yields  m^  =  -ib^. 
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In  principle,  in  order  to  solve  for  R(z),  an  additional 
measurement  would  be  required  to  obtain  the  residues  b^. 
However,  the  next  section  demonstrates  that,  in  the  presence 
cf  slow  velocity  layers  in  the  sediment,  the  frequency  can  be 
selected  low  enouqh  to  eliminate  the  trapped  modes  (or  bound 
states  of  the  Schrddinqer  equation),  and  thus  the  continuous 
part  of  the  solution  for  R(z)  will  suffice. 

4 . 4  Round  States 

The  potential  diaqram  (Fig.  3)  indicates  that  bound 
states  may  occur  due  to  the  presence  of  a  low  velocity  zone 
near  the  water  sediment  interface.  Bound  states  are  square 
inteqrable  solutions  of  Schrbdinger ' s  equation  and,  as  will 
be  shown  later,  present  considerable  difficulty  in  the 
inversion  procedure.  The  number  of  bound  states  M  was 
obtained  by  Bargmann^^ 

00 

M  <  /  z| V  (z) |dz  <  M  +  1  (4.50) 

0 

where  V_ ( z )  is  the  negative  portion  of  the  potential  for 
z  >  0.  It  is  clear  that  the  number  of  bound  states  is 
determined  by  the  width  of  the  low  velocity  zone  (prescribed 
by  the  neology)  and  by  the  depth  of  the  potential  well  which 
is  a  function  of  the  frequency  at  which  the  experiment  is 
conducted . 
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To  eliminate  bound  states,  we  impose  the  condition: 


/  z|v_(z)|dz  <  1 
0 


Now , 


/  z | V_ ( z ) | dz  < 
0 


1 

2 


V  (z) 


I  max 


(4.52) 


where  l  is  the  width  of  the  well  and  |v_(z)|  its  maximum 
depth: 


V-(z) I  max* I T" 


(4.53) 


The  condition  expressed  by  eq.  (4.12)  is  satisfied  when 


(4.54) 


In  particular,  for  — —  <  ,  there  will  be  no  trapped 

cmin 

modes.  This  simplified  condition  is  a  refinement  of 

Stickler's  result  that  (— )  should  be  sufficiently  small 

c0 

(4.50).  We  have  found  that  either  of  the  two  simplified 
conditions  is  too  restrictive  in  practice  and  one  should  use 
our  full  equation  (4.54).  Hamilton  has  studied  the  charac¬ 
teristics  of  surface  sound  channels  in  marine  sediments^®)  . 
The  velocity  ratio  R  ■  ^cmin/c0^  ranges  from  0.984  for 
pelaqic  clay  to  0.99  for  terriqeneous  sediments,  while  the 
heiqht  of  the  channel  depends  on  the  velocity  gradient  a 
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(4.55) 


H  =  c q ( 1  -  R)/a 


Therefore,  for  R  ~  1,  the  condition  of  zero  bound  states  (eq. 
4.54)  can  be  written 


u)  <  a  ( 1 


(4.56) 


A  marked  improvement  in  the  hound  obtains  if  a  linear 
velocity  orofile 


c ( z )  =  az  +  c  . 

min 


(4.57) 


with 


c(  Z ) 


c 


0 


is  assumed  in  the  condition  of  zero  bound  states  (eq.  4.51). 
The  correspondinq  potential  is  then 


V(z) 


f  az  +  c  .  ) 
min-' 


for 


0  <  z  <  i 


(4.58) 


■Substitutinq  this  V(z)  into  eq.  4.51  and  inteqratinq  by  parts 
yields 
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(4.59) 


2  1  o 

-  [ ( 1 -R )  +  j  (1-R)  +  in  R]  <1 

a 

For  R  ~  1, 

-  in  R  =  (1-R)  +  j  (1-R)2  +  -j  (1-R)3  (4.60) 

Therefore , 

u  <  a/1  ( 1-R ) ” 3/2  (4.61) 

Table  I  presents  the  upper  bound  on  the  probing 
frequency  with  the  condition  of  zero  bound  states  for  repre¬ 
sentative  values*  of  a  and  R  in  the  abyssal  plain 
environment^^.  The  current  frequency  of  operation  ^  2  ^ , 

220  Hz,  is  low  enough  to  eliminate  the  bound  states  in  clayey 
silt  and  silty  clay.  It  is  assumed  that  for  operations  in 
clay  sediments,  an  acoustic  source  will  be  available  at  about 
half  the  current  frequency  which  would  be  sufficient  to  do 
away  with  the  possible  bound  states. 


The  velocity  gradient  (a)  can  assume  values  over  a  wider 
range  than  shown  in  Table  I.  For  instance,  Frisk  et 
al.,'4'  have  inferred  from  experimental  data  that  (a) 
ranged  for  0,5s-1  to  2.9s-1  at  three  locations  in  the 
Icelandic  Basin. 
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Table  I:  Upper  bound  on  the  probing  frequency  for 
the  condition  of  zero  bound  states  as  a 
function  of  velocity  ratio,  R,  and 
velocity  gradient,  a. 


^ — 

\  R 

•\ 

Clayey  Silt 

0.999 

Silty  Clay 

0.99 

Clay 

0.984 

f 

Is-1 

8.7  kHz 

275  Hz 

136  Hz 

1.2s"1 

10.4  kHz 

330  Hz 

162  Hz 

1.3s"1 

11.3  kHz 

J - 

358  kHz 

176  Hz 
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CHAPTER  V 


NUMERICAL  SOLUTION  OF  THE  GELFAND-LEVITAN  EQUATION 

5 . 1  Introduction 
The  Gelfand-Levitan  Equation 

fZ 

R(z+y)  +  K(z,y)  +  J  K(z,z')  R(y+z')  dz ' 

—  00 

is  a  Fredholm  equation  of  the  second  kind  in  the  variable  y 
with  z  regarded  as  a  parameter.  The  object  of  this  chapter  is 
to  present  three  methods  of  computing  K(z,z)  and  hence  to 
reconstruct  the  potential  V(z). 

5 . 2  Series  Expansion 

A  parameter  A  is  introduced  in  the  Gelfand-Levitan 
equation 

rZ 

R(z+y)  +  K(z,y)  +  A  K(z,z')  R(y+z')dz'  (5.2) 

-y 

and  a  solution  is  sought  by  a  method  of  successive  approxima¬ 
tions.  The  solution  is  written  as  a  power  series  in  A 

K(z,y)  =  Kq ( z ,y)  +  AK1(z,y)  +  A2K2(z,y)  +  ... 

(5.3) 


=  0 

(5.1) 
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Integrating  term  by  term  and  equating  coefficients  of  equal 
power  of  X  one  gets 

Kq ( z , y )  =  -R (z+y)  (5.4) 

K-^Zfy)  =  -  jZ  R(y+z ' )  K0(z,z')dz'  (5.5) 

-y 

and  in  general 

Kn(z,y)  =  -  j  R(y+z ' )  Kn_1 (z , z ' ) dz '  (5.6) 

-y 

The  theory  of  Volterra  integral  equations  of  the  second 
kind  demonstrates  that  the  series  is  convergent  for  all  X 
when  the  norm  of  R, 

( ^  o 

||R| I  *  j  R  (z+y)dy  (5.7) 

-z 

( 6  i ) 

exists.  But,  Parseval ' s  identity  indicates  that  | ( R | |  is 
always  finite  since 

✓  00  00 

I  I R I  l  <  j  R2  (y)dy  =  |b(p)|2dp  (5.8) 

—  00  —CO 

2 

Hence,  b(p),  given  |b(p)j  £l  and  b(p)^>  (1/p  ),  is  square 
integrable . 
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The  Kn(z,y)  are  bounded  by 


Kn(z,y)  I  <  I  |R| 


*  (n-l)J 


(5.9) 


where  A  = 


Z  ( Z  ( z 

R2  (y+z  '  )dz  '  and  B  =  dy  'dy"R2  (y '  +y" ) 


-y 


-y  -y 


5.2.1  The  Potential 

The  expansion  of  K(z,y)  yields  a  corresponding  expansion 
for  the  potential  V(z)  =2  -K^fZ^  , 


V(z) 


=  +  V*2)  +  ... 


(5.10) 


where , 


V(0)  =  2 


dK 


dz 


(z,z)  = 


-2  §-(2z) 


(1)  ^1  2 

=  2  -^4-(z,z)  =  4RZ  (2z ) 


dz 


(2) 


(5.11) 


(5.12) 


The  computation  of  V  is  more  involved  and  is  presented 
here  for  reference  purposes 


K2(z,y)  = 


R (y+z ' )  R(z'+z")  R (z+z" ) dz"dz ' 


-y  -z  ' 


(5.13) 
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which  can  be  written 
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Finally, 


V(2)  (z)  =  4R(2z) 


2z 


R2 (z ' ) dz  '  + 


+  2 


J 

-z 


R(z+z  ' )  R(z 


+z")  — (z+z ' ) dz "dz ' 

a  Z 


-Z 


(5.20) 


V(2) (z)  displays  the  global  character  of  the  higher  order 
terms  (which  become  increasingly  unwieldy) . 

To  summarize,  the  following  approximation  will  be  used 


V  (z) 


9dR 

dz 


(2z)  +  4R2(2z) 


(5.21) 


5.2.2  Connection  with  Other  Formulations 
To  first  order, 

v(z)  =-2~(2z)  (5.22) 

which  can  be  written  in  terms  of  the  reflection  coefficient 
b  (p) 


V(z) 


2i 

Tf 


(  pb (p) e  i2pzdp. 


(5.23) 


Recall  the  expression  for  the  reflection  coefficient 
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b(p)  =  eipz  v(z')  $(z')dz' 


(5.24) 


The  Born  approximation  <J>(z')  -  eipz  results  in 


b{p)  *  —  V(z')e2ipz'dz 


(5.25) 


Taking  a  Fourier  transform  of  both  sides  leads  to 


V(z) 


=  %  pb(p)e  "2ipzdp 


(5.26) 


That  is  the  result  obtained  by  Cohen  and  Bleistein  (45) . 

The  first  order  approximation  to  V(z)  in  the  Gelf and-Levitan 
formulation  is  therefore  the  Born  approximation. 

Stickler's  inverse  procedure  entails  the  solution 
of  the  coupled  equations 


-u“  -  q(z)u2  =  p  u2 

‘  f0 

q(z)  =  j  pb  (p)  u2  (z  ,p)  dp 


(5.27) 


The  auxiliary  potential  q(z)  is  identical  to  V(z)  and 
^(Zjp)  is  the  solution  that  satisfies  ^(z,?)  =  e  ipz,  z-»— 00 
To  first  order  assume  u^tZjp)  =  e  lzp  and  obtain  once  again, 


91 


q(z)  =  ~  pb(p)e"2ipzdp 


(5.28) 


5.2.3  Algorithm 

The  first  step  is,  of  course,  the  computation  of  the 
Fourier  transform  of  the  data, 


R(z)  =  ^ 


b(p)e‘ipzdp, 


(5.29) 


followed  by  a  straightforward  computation  of 


V(Z)  (2Z)  +  4R  <2Z) 


(5.30) 


Observe  that  determining  the  potential  at  depth  z  requires 
data  at  depth  2z.  The  truncation  of  the  series  leads  to 
a  deterioration  of  the  estimate  of  V(z)  with  depth.  In 
particular,  since  b(p)  is  integrable,  R(z)  tends  to  zero 
as  z-*00  (Riemann-Lebesgue  theorem)  and  it  follows  that 
V(z)  -*-0  at  depth  even  when  the  potential  V(z)  tends  to  a  non¬ 
zero  final  potential  V^. 
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5 . 3  Uniqueness  of  the  Solution 

Suppose  there  were  two  solutions  of  the  Gelfand-Levitan 
equation,  K(z,y)  and  K'(z,y) 


K(z,y)  +  R(z+y)  +  j  K(z,zf)  R(y  +z ' )  dz ' 

-y 

z 

K ' (z ,y )  +  R(z+y)  +  K'(z,z')  R(ytz')  dz 

-y 


=  o 

=  o 

(5.31) 


Then  K(z,y)  =  K(z,y)  -  K' (z,y)  satisfies  a  homogeneous 
Volterra  equation  of  the  second  kind 

Z  /v 

K (z,y)  =  K(z,z')  R(y+z')  dz'  (5.32) 

-y 

Since  R(z)  is  square  integrable (cf .  Ch.III),  it  follows  (61) 

/s 

that  this  equation  has  only  the  trivial  solution  K(z,y)  =  0 
and  therefore  the  solution  of  the  Gelfand-Levitan  equation 
is  unique. 

5 . 4  Finite  Difference  Methods 

The  natural  way  to  solve  the  Gelfand-Levitan 
equation  is  to  discretize  it  by  converting  the  integral  into 
a  sum 


n 

K (n,m)  +  R(n+m)  +  h  l  w,K(n,i)  R(m+i) 

i=-m  1 


0 

(5.33) 
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where  w^  are  the  weights  in  an  appropriate  quadrature 
formula.  For  the  trapezoidal  rule,  w^  is  \  at  the  end  points 
and  1  in  between 

n 

K(n,m)  +  R(n+m)  *■  h  E  K(n,i)  R(i+m) 

i=-m 

-  ^h (R ( 0 ) K  +  R  (n+m)  K(n,n))  =  0  (5.34) 

and  since  R(0)  =  0 

n  h 

K(n,m)  +  h  E  K(n,i)  R(i+m)  +  R(n+m)(l  -  ^(n,n)]=0 
i=l-m  1 

(5.35) 

To  solve  for  K(n,n),  one  has  to  solve  the  Ware  and  Aki  (21  ) 
type  matrix  equation, 

(I  +  hR) K  =  h  (5.36) 

A 

(I  +  hR)  can  be  inverted  by  Gauss  elimination.  A  more 
efficient  algorithm  has  been  presented  by  Berryman  and 
Greene  (  26)  who  put  the  equation  in  the  form: 
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I  +  h 


Rl~ 

R1R2 

• 

hk  (n,  -n+1) 

0 

• 

R1R2*  *  *  R2n_ 

hk (n,n-l) 

1+hk (n,n) 

*  H*  O  • 

t 

(5.37) 

where  k(n,m)  *  Hl^m) -  (5.38) 

l-|K(n,n) 

To  obtain  the  potential  at  depth  n  -  *5, 

qn_h  *  ^[K(n,n)  -  K(n-1,  n-1) ]  (5.39) 

one  needs  to  invert  a  (2n-2)  x  (2n-2)  matrix  for  K(n-l,n-l) 
followed  by  the  inversion  of  a  (2n)  x  (2n)  matrix  for 
K(n,n).  In  fact,  since  the  object  is  to  reconstruct  the 
potential  down  to  depth  nh,  a  succession  of  matrices  of 
increasing  size  have  to  be  inverted.  The  Berryman  and  Greene 
algorithm  is  similar  to  the  Levinson  algorithm  for  the 
inversion  of  Toeplitz  matrices  (22  ) .  The  method  proceeds 
by  recursion;  given  the  solution  of  the  (2n-l)  x  (2n-l) 
system,  the  solution  of  the  (2n)  x  (2n)  system  is  generated 
by  using  the  recursion  formulae: 
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fi(2n) 


gi(2n) 


1  r. 


rn  1 


f i (2n-l) 


g± (2n-l) 


(5.40) 


The  vectors  f(2n-l)  and  g(2n-l)  are  defined  by 


R1R2 


f1(2n-l) 

R1 ' * ' R2n-lJ  |_f0(2n-l) 


f2„.i(2n-l) 


f2„.i(2n-l) 


g-L  (2n-l) 
gQ (2n-l) 


(5.41) 


R1R2 


Rl-  "  R2n-1 


g2n-l(2n-1)"l  rf2n-l(2n-1> 


gn(2n-!) 


fL(2n-l) 


f0(2"-1»-°2n-a 


where  Q2n  =  Q2n.1U-rn  ) 


(5.42) 


2n-l 

and  r-  =  h  Z  f . (2n-l)  R 
2n  i=0  1 


2n+l-i/Q2n-l 


The  recursion  starts  with  fQ(0)  =  1  and  gQ(0)  =  AR^ 
The  solution  is  then, 
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K (n , n) 


(5.43) 


2  1  ~  Q2n-1 (1  +  r2n) 
1  +  Q2n-l(1  +  r2n) 


Berryman  and  Greene  have  suggested  smoothing  R(z)  to  obtain 
R(n)  (26). 


R(n)  =  yi  R(z)dz 


(5.44) 


which  would  allow  for  delta  functions  in  R(z) . 


5 . 5  Coupled  Equations  Method 

A  method  to  bypass  the  matrix  inversion  has  been 
suggested  by  Kritikos,  Jaggard  and  Ge  who  were  interested 
in  determining  the  dielectric  permittivity  of  a  slab  from 
reflection  measurements  at  normal  incidence.  The  algorithm 
was  tested  numerically  on  reflection  coefficients  that 
could  be  represented  by  two  and  three  pole  Butterworth 
filters.  The  scheme  uses  in  conjunction  with  the  Gelfand- 
Levitan  equation,  the  hyperbolic  partial  differential 
equation  satisfied  by  K(z,y)  (cf.  Eq.  (4.44)). 


2 

,  v  8  K  ,  ,  nv/  \  dK(z,z) 
(z,y)  -  -2  (z,y)  -  2K (z,y)  — ~^~L 

3y 


(5.45) 
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For  simplicity,  a  change  or  coordinates  is  introduced 


and  the  equations  are  discretized: 


(5.46) 


m-1 

Gelfand-Levitan :  K  +  2h  I  K. 


m,n 


i=2 


i+n-l,m-i+l  2i-l 


+  (1  +  hKm+n-l,l)R2m-l  =  0 


(5.47) 


PDE:  K.  ..  =  K  , «  +  K  . ,  +  2h [  (K  .  ,  -  K  ,)-l]K 

m+l,n+l  m,n+l  m+l,n  m+1,1  m,l  1  m,n 


(5.48) 


The  potential  is  obtained  through 

V  ,/>==  \  (K  .  -  K  .  .  ) 
m-V^  h  m,  1  m-1, 1 


(5.49) 


The  point  of  the  method  is  that  Km  ^  can  be  computed 
directly  from  the  Gelfand-Levitan  equation  without  a 
matrix  inversion  since  the  terms  within  the  sum  can  be 
generated  via  the  PDE. 

5.5.1  Analysis  of  Stability 

Although  it  has  not  been  possible  to  study  the 
stability  of  the  coupled  system  of  equations,  a  simplified 
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analysis  of  the  PDE  may  be  of  interest. 

Let  K  =  W™ elkn  be  a  Fourier  component  and 
m,  n 

assume  that 


2h(K 


m+1 , 1 


-  K 


m.  1 


)  = 


h2V 


m 


is  known.  The  PDE  becomes 


(5.50) 


w”+1  eik  =  wVk  +  w1"*1  ♦  (h2V  -1) 

m 


(5.51) 


The  amplification  factor  is  therefore 


h2V_ 


1  + 


m 


eik-l 


(5.52) 


which  is  larger  than  1  for  arbitrarily  small  step  size  h  or 
potential  as  k  ->■  0. 

5 . 6  Error  Estimates 
5.6.1  Discretization  Errors 

The  discretization  of  the  Gelfand  Levitan  equation 
is  accompanied  by  an  error  ^  , 


n 

K(n,m)  +  R(n+m)  +  h  2  w.K(n,i)  R(m+i)  +  e  =  0 

■  l  n  f  m 

i=-m 


(5.53) 
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For  the  trapezoidal  rule,  c  =  -jj  (KR) 


The  equation  that  is  actually  solved  is 


K(n,m)  +  R(n+m)  +  h  E  w.K(n,i)  R(m+i) 

i=-m  1 


Or,  in  matrix  notation  (see  equation  5.37  ). 


(5.54) 


(I  +  hR)  K  =  b 


(5.55) 


rather  than 


(I  +  hR)K  =  b  -  e 


(5.56) 


It  follows  that 


(I  +  hR)  (K  -  K)  =  e 


(5.57) 


For  a  nonsingular  matrix  (I  +  hR) , 


K  -  K  =  (I  +  hR)  e 


(5.58) 


Taking  norms. 


K  -  K|  I  <  I  |  (I  +  hR)  ■L  .  e 


(5.59) 
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Assuming  that  J  | hR| |  <  1,  where  the  norm  is  such  that 

II  IM  =1, 


M  (I  +  hR)  -1 1  |  <_  —  -|-|  |  (5.60) 

Therefore, 

ii*  -K'i  <  i-h n r| i  <5-61» 

2n 

Using  the  «>  norm,  | | e | |  *  max  |e.|  and  ||r||  =  Z  |R(i)| 

i-1 

(which  also  equals  ||R||  ,  the  absolute  column  sum  of  column 
(2n) )  . 


1  -  h  z  | R(i) | 
i=l 


(5.62) 


2 

Since  |e|  is  0 (h  )  for  the  trapezoidal  rule,  the  approximate 

A 

solution  K  converges  to  the  exact  solution  K  as  h  -*•  0 .  Note 
that  the  denominator  in  (equation  5.62)  decreases  with 

A 

increasing  depth,  hence  raising  the  bound  on  | |K  -  k| | . 


5.6.2  Data  Errors 

Errors  in  the  plane  wave  reflection  coefficient 
r(kz)  and  in  its  Fourier  transform  R(z) ,  mean  that  one  is 
solving  for  K(n,m)  in  the  Gelf and-Levitan  equation: 
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K (n,m)  +  R (n+m)  +  h  Z  w.  K(n,i)  R(m+i)  =  0 

i=-m  1 


(5.63) 


Or  in  matrix  notation. 


A  A  A 

(I  +  hR)  K  =  b 


which  can  be  written. 


(I  +  hR  +  h6R)  (K  +  6  K)  =  (b  +  6b) 


(5.64) 


(5.65) 


A  A 

where  R  =  R  +  6R,  and  K  =  K  +  6K. 
It  follows  that 


6K  =  (I  +  (I  +  hR)_1  h6R) _1  (I  +  hR)  1  (6b  -  Kh6R) 


Hence, 


(5.66) 


| ] 6K||  £  | | (I  +  hr)"1 j  | . | | (I  +  (i  +  hR)"1  h6R)-1| | 

•  ( I  I  ^b |  |  +  h||6Rj  |  .  |  |K|  | ) 

(5.67) 

or, 

ll6Kll  1  r-  h(|  |R|  I'  '+  |  1  6R|  [  )  (M6bM  +  h|  |  K  |  |  .1  |  6R|  I  ) 

(5.68) 
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As  the  errors  in  R(z)  increase,  the  bound  on  ||6k|| 

also  increases.  Moreover,  it  is  clear  that  as  depth 

2n 

increases  =  £  |R(i)|  also  increases  raising  the 

i=l 

bound  on  the  relative  error  | | 6K | | / | | K | | . 

5.6.3  Errors  in  the  Potential 

The  discrete  version  of  the  potential  yields: 

j  6 V ( n )  i  <£  ( | 6K (n , n) |  +  |6K(n~l,  n-1) | )  (5.69) 

In  the  absence  of  data  error  6R,  )6V(n)|  ->  0  as  the 

sampling  interval  h  -*■  0  since  the  quadrature  error  e  is 
2 

0 (h  )  (see  equation  5-62). 
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CHAPTER  VI 


NIJM E RICAL  RESULTS 


The  derivation  of  the  plane-wave  reflection  coefficient 
from  pressure  measurements  for  the  setup  of  the  Frisk  experi¬ 
ment  was  the  subject  of  Mook's  thesis^^.  Unfortunately, 
the  attempt  by  Mook  to  apply  his  technique  to  one  set  of 
experimental  data  yielded  reflection  coefficients  higher  than 
one;  such  a  result  is  not  consiscent  with  physical  require¬ 
ments  ( | P (  <  1).  It  is  not  yet  clear  whether  Mook’s  problem 
lay  with  improper  modeling  of  the  experimental  setup,  with 
the  imprecision  of  the  data,  or  with  the  numerical  techniques 
used  to  extract  the  reflection  coefficient  from  the  data.  It 
should  be  noted,  however,  that  the  experimental  data  have 
yielded,  via  trial  and  error  methods,  excellent  models  for 
the  acoustic  parameters  of  the  seabed^^.  This  was  done  by 
assuming  a  seabed  model,  comouting  the  pressure  field  and 
matching  it  up  to  the  observations. 

The  starting  point  in  our  analysis  is  a  layered  model  of 
the  acoustic  profile  to  be  recovered.  From  this  exactly 
known  model,  a  plane-wave  reflection  coefficient  is 
generated.  The  reflection  coefficient  is  what  would  have 
been  computed  from  Mook's  method,  or  subsequent  improvements 
to  it,  for  the  Frisk  experiment  on  this  particular  layered 
seabed.  The  reflection  coefficient  is  then  Fourier  trans- 
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formed  via  an  FFT  to  provide  the  input  to  the  Gelfand-Levitan 
integral  equation  in  accordance  with  the  theory  developed  in 
Chapter  IV.  The  particular  case  of  two  half-spaces  in 
contact  (step  discontinuity)  yields  a  reflection  coefficient 
which  can  be  Fourier  transformed  analytically.  Closed-form 
solutions  for  the  Born  and  Improved  Born  approximations  to 
the  inverse  solution  are  also  obtained  for  that  case. 

More  generally,  the  Gelfand-Levitan  equation  is  solved 
numerically  to  yield  the  potential  and  velocity  profiles. 

We  have  concentrated  our  efforts  on  the  study  of  the 
Gelfand-Levitan  inverse  method  using  synthetic  data  for  which 
the  correct  answer  is  known.  The  comparison  of  the  recon¬ 
structed  profile  with  the  known  original  profile  allows  us  to 
assess  the  impacts  of  limited  aperture,  frequency,  profile, 
noise,  density  and  path  loss  on  the  accuracy  of  the  numerical 
schemes  described  in  Chapter  V, 

6 . 1  Generation  of  the  Reflection  Coefficient 

The  first  step  in  the  evaluation  of  our  approach  is  the 
generation  of  the  reflection  coefficient.  With  the  exception 
of  the  simplest  cases,  the  plane-wave  reflection  coefficient 
must  be  generated  numerically.  The  method  used  is  based  on 
the  Thomson-Haskel 1  propagator  matrix  approach  used  by 
Mook^'U,  with  one  modification.  The  plane-wave  reflection 
coefficient  is  obtained  as  a  function  of  the  vertical  wave- 
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number  (R(kz))  rather  than  as  a  function  of  the  horizontal 
wavenumber  (R(kx)).  There  are  two  essential  elements  of  the 
propagator  matrix  approach:  (a)  dividing  the  acoustic 
profile  into  homogeneous  layers,  and  (b)  selecting  variables 
that  are  continuous  across  interfaces.  The  latter  is  accom¬ 
plished  by  choosing  the  pressure  P ( z  )  and  vertical  component 
of  velocity  U(z).  Within  an  isovelocity  layer,  the  field  can 
be  decomposed  into  up  and  down  going  waves, 

ik  z  -ik  z 

P(z )  =  P+  e  2  +  P_  e  2  (6.1) 


where  kz  denotes  the  vertical  wavenumber  within  the  homogene¬ 
ous  layer.  Unlike  P(z),  P+  and  P_  are  discontinuous  at  an 
interface.  The  normal  component  of  velocity  U(z)  is  related 
to  P(z)  through  one  of  the  time-harmonic  "telegraph" 
equations : 


3P 

3z 


iaipu 


(6.2) 


which  yields. 


U(z) 


ikzz 


P 


-ik 

e 


(6.3) 


Where  Y  =  k  /up  is  the  admittance  of  the  homoqeneous 
n  z  n  3 

n 

layer  n.  In  matrix  form,  equations  (6.2)  and  (6.3)  become: 
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P(z  ) 

U(Z  ) 


ik  z  -  ik  z 

A  z  z 

6  e 

P^ 

+ 

ik  z  -ik  z 

Y  e  z  -Ye  z 

P_ 

n  n 

. 

(6.4) 


When  P(z)  and  U(z)  are  known  at  point  z\  of  the  layer, 
one  can  deduce  the  up  and  down  going  components  P+  and  P_  by 
inverting  the  matrix  in  (6.4) 


P+ 

r.ik^i 

e 

-lkzzl  1 
e 

-1 

P_ 

'lkzzl 

|_-v 

“ikzzl 
-Ye  z  1 

n  _ 

(Within  a  homogeneous  layer  P+  and  P_  are,  of  course, 
constants.)  One  can  then  determine  P(z)  and  U(z)  at  another 
point,  say  z2,  within  the  layer  by  substituting  (6.5)  into 
the  right-hand  side  of  (6.4): 


p(zl) 

coskz(z1  -  z2) 

— 

Y"  sinkz(zx  -  z2) 
n 

p(z2) 

U(Z1J 

iY  sinkz(zx  -  z2) 

coskz(z1  -  z2) 

U(z2) 

(6.6) 
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or,  introducing  the  propagator  matrix  $ 


p(*l) 

p( 2  2) 

»— -i 

N 

1 

CM 

N 

& 

II 

U(Zl) 

U(z2) 

One  can  now  proceed  from  one  layer  to  another  by  integration, 


0 

N 

lE _ 

=  v*n 

_P(2n) 

Lu(zo) 

* 

U(zn)_ 

Where  the  depend  on  the  parameters  of  the  material 

making  up  the  respective  layers.  Mook^64^  has  found  that  the 
computations  can  be  improved  by  modifying  (6.8)  to 


_p(zo) 

n 

=  n 

_V(zo)_ 

i=l 

(6.9) 


in  which 


n 

n 

" a  i 

bi 

11 

♦ll 

♦l2 

i  =  l 

b‘ 

ci  ai_ 

__*21 

*22 

(6.9a) 


-108- 


and  where 


^  n  1 

=  — y —  is  a  normalized  admittance 


a.  =  cos  k  h. 
i  z.  i 


b.  =  -l  sin  k  h. 
l  z  .  1 

l 


h  • 


th 


^  =  thickness  of  i  layer 


kZl  *  'kn  '  (k0  '  kZ()) 


This  has  the  effect  of  giving  the  two  components  P(z)  and 
YU(z)  similar  scales  of  magnitude. 

Now,  the  reflection  coefficient  R(kz)  at  the  top  inter¬ 
face  is  defined  as 


Therefore  by  introducing  (6.5)  into  (6.9)  one  obtains 


1 

1 

1 

*11 

♦l2 

p(*n) 

2 

1 

-1 

♦21 

*22 

'nU(2n) 

_  _ 

(6.11) 
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To  compute  the  ratio  in  (6.10),  a  relation  (boundary  condi¬ 
tion)  is  needed  between  P(zn)  and  U(zn).  That  relation  is 
obtained  by  noticing  that  in  the  (n  +  l)st  layer,  one  has 
only  down  going  waves.  Therefore, 


Kfk  1  =  *n~  +21*  cN+l[*12~  *22^ 

Z  P+0  *11+"  *21  +  ?N+1^12+  *22^ 


(6.14) 


Mook^^  has  also  found  it  advantageous  to  scale  the 
layer  propagation  matrices  so  that  the  largest  element  value 
in  a  given  layer  matrix  is  1. 

The  values  of  admittance  for  three  terminations  are 
readily  identified:  hard  bottom  (Yn  +  ^  =  0),  soft  bottom 
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(Yn+1  =  » ) ,  and  isovelocity  termination  (Yn+1  =  Yn ) . 

Although  of  theoretical  interest,  and  of  possible  applicabil¬ 
ity  to  other  situations,  the  soft  bottom  case  has  no  appli¬ 
cation  in  the  ocean  bottom  problem  discussed  here. 

The  division  of  the  acoustic  profile  into  homogeneous 
layers  is  done  in  a  way  akin  to  quadrature  formulae  in 
numerical  integration  using  thinner,  closely-spaced,  layers 
in  regions  of  rapid  change  in  the  acoustic  parameters  and 
thicker,  wider-spaced  layers  in  regions  of  slow  change.  We 
have  found  that  the  common  rule  of  thumb,  ten  layers  per 
wavelength,  although  satisfactory  in  general,  is  probably  too 
conservative.  We  have  found  that  for  complicated  profiles, 
where  the  reflection  coefficient  was  required  for  a  large 
number  of  values  of  kz ,  the  computation  of  the  plane-wave 
reflection  coefficient  constituted  the  most  time-consuming 
step  in  modeling  the  whole  inversion  procedure.  Clearly, 
this  step,  is  inherent  only  in  our  analytical  evaluation  of 
the  inversion  problem;  when  reflection  coefficients  are  being 
processed  from  measurements,  this  step  will  be  eliminated. 

In  order  to  test  our  inverse  procedure,  we  selected  a 
few  representative  profiles  which  are  defined  as  follows: 
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(a)  Half-Space 


The  half-space  profile  is  defined  by  a  velocity  profile 

c(z  ) , 


c (z )  =  cQ  , 

=  C]  , 

The  reflection  coefficient  for  this  profile  and  its 
Fourier  transform  were  obtained  analytically  in  Section 
(6.2). 

(b)  Window  Profile 

Although  of  no  direct  application  in  ocean  bottom 
acoustics,  the  window  profile  defined  by 

c(z )  =  c0  ,  z  <  0 

=  Cl  ,  0  <  z  <  L 

=  c0  ,  z  >  L 

is  a  standard  example  used  to  test  inverse  procedures.  The 
window  profile  is  of  particular  interest  in  electromagnetics 
(dielectric  slabs). 


z  <  0 
z  >  0 
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(c)  Stickler  Profile 


Stickler^"^  chose  a  twice  continuous  function  to  test 
the  Deif t-Trubowitz  inverse  method, 

z  <  0 
0  <  z  <  L 
L  <  z 

This  profile  does  not  support  a  trapped  mode  at  any 
frequency. 

Stickler  generated  the  associated  reflection  coefficient 
via  the  Riccati  equation,  while  we  used  the  Thomas-Haskel 1 
procedure  outlined  earlier  in  this  section.  An  example  of 
the  Stickler  profile  is  shown  in  Fig.  10. 

(d)  Frisk  Profile 

This  velocity  profile  is  based  on  the  results  obtained 
by  Frisk  in  his  deep  ocean  bottom  experiments^ 

c(z )  =  c0  ,  z  <  0 

=  Cq  -  (co-ci )z  ,  0  <  z  z  1 

=  Cj  +  ,  1  <  z  <  L 

=  c  2  ’  z  >  L 


=  cQ  +  (crc0)(3(f-)  -2(g)  ), 

=  Ci 
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Typical  parameters  used  are:  cq  =  1540  m/s,  = 

1515  m/s,  C2  =  1655  m/s,  y  =  .97,  and  L  =  145  m  (see 
Fiq.  11). 

In  spite  of  the  low  velocity  zone  near  the  ocean-bottom 
interface,  this  profile  does  not  support  trapped  modes  at  the 
frequency  used  for  the  examples  (25  Hz). 

Examples  of  reflection  coefficients  are  shown  in  the 
amplitude  and  phase  diagrams  of  Fig.  18.  Note  that  as  we 
have  shown  in  Chapter  III  for  the  uniform  density  case,  the 
amplitude  goes  asymptotically  to  zero  as  the  vertical  wave- 
number  goes  to  infinity. 

6 • 2  Case  of  a  Step  Discontinuity  in  Potential 

We  have  been  able  to  obtain  analytically  the  Fourier 
transform  of  the  reflection  coefficient  in  the  case  of  a  step 
discontinuity  in  potential,  that  is,  in  the  case  where  the 
ocean  bottom  is  a  homogeneous  half-space.  (V(z)  =  Vq  for 
z  >  0) . 

The  reflection  coefficient  at  a  step  discontinuity  in 
potential  is 
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r(k7.)  = 


kz~  ^~V0  +  kz 
kz  *  kz: 


[k,-  /-vn*  k^l2 


V, 


The  Fourier  transform  of  r(k„)  is 


(6.15) 


R(z)  = 


77 


-ik  z 

r(k  )e  dk 


(6.16) 


i«  (s  -  /Vq+  s  ) 2 
-il  VQ 


sz 


ds 


which  in  this  form  can  be  identified  as  a  known  inverse 
Laplace  transf orm ^ ^ ) , 

R(z)  =  -  |  J2  (/VJ*)  (6.17) 

Thus,  in  this  particular  case,  one  can  proceed  to  the 
Gelf and-Levitan  procedure  with  an  input  which  is  as  accurate 
as  the  computation  of  the  Bessel  function,  and  therefore, 
this  approach  removes  any  inaccuracies  which  may  be 
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introduced  either  in  the  computation  of  the  reflection 
coefficient  or  in  that  of  its  Fourier  transform. 

Figs.  12a, b,c  represent  the  reflection  coefficient  of 
eq.  6.15  and  its  Fourier  transform  for  a  half-space 
( eg .  6.17). 

As  a  matter  of  fact,  we  have  illustrated  in  Figs.  13  and 
14b  reconstructions  of  the  step  discontinuity  using  the  exact 
analytical  expression  (eq.  6.17)  for  the  Fourier  transform  of 
the  reflection  coefficient  and  the  numerically  devised 
transform  in  accordance  with  the  method  discussed  in  Section 
6.1.  It  can  be  seen  that  our  reflection  coefficient  method 
yields  results  as  accurate  as  the  Bessel  function  expression. 

6.2.1  Approximate  Solution  of  the  Gelf and-Levitan  Equation 
for  a  Step  Discontinuity 

Continuing  with  the  case  of  a  homogenous  half-space,  or 
step  discontinuity,  the  "Improved  Born"  approximation 
presented  in  Chapter  V  can  be  now  readily  obtained  in  terms 
of  Bessel  functions 
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V(z  ) 


=  -  (2z)  +  4R2(2z) 


(6.18) 


✓vn 

-  4—  ^C^V)  -  h  j2(2/V)  -  T-  j2(2/V) 

z  z 


Born 


Correction  Term 


Improved  Born 


From  the  asymptotic  behavior  of  the  Bessel  functions,  we  know 
that. 


(6.19) 


Hence,  we  can  deduce  the  asymptotic  behavior  of  the  recon¬ 
structed  potential, 


V(z)  -►  0  as(z’3^2) 

z  +  * 


(6.20) 


which  recalls  the  limitation  of  the  Improved 
mation,  already  mentioned  in  Chapter  V,  that 
in  the  case  of  a  finite  terminal  potential. 


Born  approxi- 

A 

V(z )  0,  even 

Z  « 
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On  the  other  hand,  the  asymptotic  behavior  of  the  Bessel 


functions  near  the  oriqin  is 


Ji  U )  ♦ 

c  -  0 


1 

7  ? 


(6.21) 


J?  (?  )  + 

c  -  0 


1  2 
7  5 


and  therefore,  substituting  in  (6.18)  the  reconstructed 

A 

potential  v(c  )  tends  to  the  exact  potential  near  the  origin, 


V( z  )  -  V  (6.22) 

z  -*■  0  u 

At  the  origin,  the  "Improved  Born"  correction  term  is  zero, 
but  its  contribution  to  the  accuracy  of  the  results  becomes 
progressively  more  important  as  z  increases. 

We  illustrate  the  Improved  Born  approximation  in 
Fig.  12d  along  with  the  Standard  Born  approximation  and  the 
Correction  Term.  The  substantial  improvement  due  to  the 
Correction  Term  in  (6.18)  over  the  standard  Born  approxi¬ 
mation  is  clearly  visible  in  the  half-space  case.  Moreover, 
our  Improved  Born  approximation  results  in  a  more  accurate 
reconstruction  of  the  acoustic  velocity  profile  to  further 
depth.  The  other  features  of  the  approximation  are  also 
visible,  i.e.,  excellent  reconstruction  near  the  origin  and 
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deterioration  of  the  reconstructed  velocity  with  depth  (the 
potential  V  4-  0  as  z  -►  »  means  that  c(z)  -►  cq,  the 

water  acoustic  velocity  as  z  ->•  ®).  The  Improved  Born  is 

also  illustrated  in  the  case  of  a  Stickler  profile  in  Figs. 

28  and  29. 

6 . 3  Fourier  Transform  of  the  Reflection  Coefficient 

The  Gelf and-Levitan  algorithm  requires  as  an  input  the 
Fourier  transform,  R(z),  of  the  reflection  coefficient  r(kz), 

R  ( z  )  =  i-  J  r(  k  )e  ^^dk  (6.23) 

2ir  1  v  z J  z 

— oo 

which  requires  a  knowledge  of  r(kz)  over  the  whole  line 
-»  <  kz  <  But  the  symmetry  property  demonstrated  in 

Chapter  III  (eq.  3.6), 

r(-kz)  -  r*(kz)  (6.24) 

reduces  the  requirement  to  a  knowlege  of  r(kz)  over  the  half¬ 
line  0  <  kz  <  “•  It  follows  that 

1  "  _ikzz 

R(z)  =  3—  /  r(k  )e  dk  +  complex  conjugate  (6.25) 

Zn  q  z  z 

is  real,  and  therefore  all  the  quantities  involved  in  the 
Gelfand-Levitan  algorithm  are  also  real. 
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Theoretically,  r(kz)  can  be  calculated  for  a  given 
acoustic  model  over  the  range  0  <  kz  <  » .  In  practice, 
samples  of  r(kz)  are  given  over  a  finite  range  a  <  kz  <  b, 
and  presumably  this  range  is  restricted  to  real  angles 
0  <  kz  <  kg,  where  kz  =  0  corresponds  to  grazing  incidence 
and  kz  =  kg  (the  water  wavenumber)  corresponds  to  vertical 
incidence.  It  is  useful  to  study  the  behavior  of  the  Fourier 
transform  under  different  restrictions  imposed  upon  the 
knowledge  of  the  reflection  coefficient  such  as  limited 
angular  aperture  and  different  sampling  densities. 

6.3.1  Fast  Fourier  Transform 

The  computation  of  the  FFT  for  reflection  coefficients 
r(kz)  computed  over  0  <  kz  <  a,  including  the  case  where 
a  >  k0  (corresponding  to  complex  angles  of  incidence),  does 
not  present  difficulties.  The  adopted  algorithm  uses  time 

decomposition  with  input  bit  reversal ( ^6  )  ^  in  fact,  for  the 

acoustic  profiles  tested,  a  =  2  or  3  kg  was  sufficient  as  the 
asymptotic  decay  of  the  reflection  coefficient  with 
increasing  kz  was  even  more  rapid  than  the  (l/kz  )  derived  in 
Chapter  III.  The  errors  in  the  imaginary  part  of  the  FFT, 
which  are  of  the  order  of  10"^  (single  precision),  are 
negligible,  and  therefore  the  real  part  of  the  FFT  and  its 
amplitude  are  interchangeable  on  the  plotted  results. 
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6.3.2  Limited  Angular  Aperture 

Since  it  is  not  practical  to  have  the  sampling  process 
cover  the  entire  non-zero  portion  of  r(kz),  it  is  useful  to 
study  the  effect  of  limited  angular  aperture,  i.e.,  r(kz) 
given  over  the  finite  range  a  <  kz  <  b,  on  R(z),  the  Fourier 
transform  of  r(kz).  The  Fourier  transform  is  affected  by  the 
nature  of  r(k_)  and  by  the  degree  of  truncation. 

For  N  sampling  points  spaced  Ak_  apart,  the  total 

sampling  interval  T  covered  is  N»Ak_.  This  corresponds  to  an 

-1  NAkz 

angular  aperture  of  a  =  sin  1  ( — r — ),  when  kz  =  0  is  within 

K0 

the  known  aperture  (i.e.,  90°  >  0  >  cos“* (Na kz /kg ) ) .  Note 
that  for  T  >  kg,  cos(T/kg)>l,  and  the  aperture  includes  all 
real  angles  plus  complex  angles.  Now,  the  Fourier  transform 
of  r(kz)  can  be  written 

.  T  -ik  z  ,  ®  -ik  z 

R(z)  =  f  r(k  le  dk  +  -n—  /  r(  k  )  e  Z  dk 

2ir  ^  z 1  z  2tt  ip  K  7,  ’  z 


=  R1 (z )  +  R2 (z ) 


(6.26) 


^l(z)  is  the  part  of  R(z)  that  is  approximated  by  the 
discrete  Fourier  transform  (N  samples). 

P'2(z)  represents  the  error  incurred  due  to  the  limited 
angular  aperture.  The  energy  lost  in  the  process. 
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(6.27) 


F  =  f  |  R-  (z  )  |  2dz 
T 

can  be  computed  on  synthetic  examples.  Due  to  the  precipi¬ 
tous  droD  in  r(kz)  bevond  kc,  the  wavenumber  corresponding  to 

—  1  c  0 

the  critical  angle  Qc=  sin  (— )  ,  and  due  to  the  asymptotic 
behavior  of  r(kz)  for  large  kz  ,  it  is  readily  shown  on 
computer  simulations  that  a  knowledge  of  r(kz)  over  real 
angles  is  adequate  for  most  cases.  In  fact,  there  is  little 
change  in  the  reconstructed  profile  as  more  angles  are 
included  beyond  the  real  ones.  On  the  other  hand,  as  the 
angular  aperture  is  restricted,  the  profile  reconstruction, 
via  the  Gelfand-Levitan  method,  produces  a  smoothed  out 
version  of  the  original  profile.  This  low-pass  filtering 
phenomenon  is  best  understood  by  interpreting  the  effect  of  a 
finite  aperture  as  a  low-pass  filtering  of  the  original 
velocity  profile.  This  may  be  seen  from  the  Rorn  approxi¬ 
mation  eq .  (5.26), 

A  .»  -2 ik  z 

V(z)  =  IT  /  kzr(k*)  e  Z  dkz 

-09 

A 

where  V(z)  and  kzr(kz)  form  a  Fourier  transform  pair; 

A 

windowing  r(kz)  signifies  low-pass  filtering  V(z)  and 
therefore  yields  a  smoothed  out  velocity  profile.  The 
Gelfand-Levitan  reconstruction  can  then  be  interpreted  as  an 
unfiltered,  faithful  reconstruction  of  the  smoothed  out 
original.  Severe  degradation  of  the  reconstruction  is 
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observed  when  the  angle  of  incidence  is  restricted  to 

90°  >  0  >  0C  (angular  aoerture  a  <  (90°-©c)).  This  is  easily 

understood,  as  the  critical  angle  region  contributes 

substantially  to  the  reflection  coefficient.  A  cursory 

examination  of  reflection  profiles  will  demonstrate  this 

point. 


A  series  of  figures  displays  the  progressive  deteriora¬ 
tion  in  profile  reconstruction  as  the  angular  range  is 


restricted.  The  evolution  from  Fig.  19  (k, 
plex  angle  of  incidence)  to  Fig.  20  (kz 


=  .512,  corn- 

max 

=  0.128),  and  Fig. 


21  (k  =  0.064,  0  =  51.13°)  where  k_  >  k  .  .  = 

max  zmax  ^critical 

0.0373  (ec  =  68.5°)  shows  that  only  small  changes  take  place 

and  that  these  changes  are  confined  for  the  most  part  to  the 

velocity  drop  region  near  the  ocean  bottom  interface. 

However,  as  the  angles  are  further  restricted  to  beyond  the 

critical  angle  (k„  <  k7  .)  major  changes  do  occur  as 

zmax  zcritical 

seen  in  Fig.  22  (k  =  .032,  0  =  71°)  and  Fig.  23  (k- 

4raax  ‘max 


0.016,  0  =  80.9°).  The  examination  of  the  impact  of  limited 


angular  aperture  in  this  case  leads  us  to  expect,  in 
practice,  a  good  reconstruction  of  acoustic  velocity  profiles 
from  reflection  data  restricted  to  real  angles  (0  <  0  <  90°). 
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6.3.3  Samplin 


The  sampling  rate  to  be  chosen  is  governed  by  three 
competing  considerations: 

(a)  Adequate  sampling  in  the  kz-domain  to  avoid  aliasing  in 
the  Fourier  transform  of  r(kz). 

(b)  Adequate  sampling  in  the  z-domain  to  obtain  a  stable  and 
accurate  Gelfand-Levi tan  numerical  reconstruction. 

(c)  Adequate  sampling  in  the  z-domain  to  obtain  the 
necessary  resolution  in  the  reconstructed  profile. 

Each  of  these  points  is  discussed  next. 

(a)  Aliasing 

Since  the  Fourier  transform  of  the  reflection  coeffi¬ 
cient,  as  is  the  rule  with  spectra  of  transients,  tends  to  be 
smooth  and  approaches  zero  asymptotically  as  kz  increases  to 
infinity,  the  sampling  interval  Akz  is  chosen  so  that 
essentially  all,  rather  than  all,  the  spectral  content  of  the 
waveform  is  contained  below  l/(2Akz).  For  a  box-like 
reflection  coefficient  of  width  2»kcr^t^ca^  where  ^critical 
corresponds  to  0crit'  the  bandwidth  is  proportional  to 
1/(2  kcritical)  which  indicates  that  an  appropriate  sampling 
interval  should  be  a  fraction  of  kcr£tical*  0ne  should  note 
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however,  that  in  the  context  of  an  actual  experiment,  the 
choice  of  Akz  is  not  as  straightforward.  The  reflection 
coefficient  r(kz)  is  obtained  from  r(kr),  the  reflection 
coefficient  as  a  function  of  horizontal  wavenumber  kr, 


(6.28) 


The  sampling  interval  in  kr  is  determined  by  the  maximum 
distance,  D,  between  source  and  receiver  during  the 
experiment 


Akr  =  i  .  (6.29) 

It  is  clear  that  a  uniform  sampling  in  kr  does  not  lead 
to  a  uniform  sampling  in  kz , 

*  -  tan  ainc  '  (6'30) 

r 

which  also  shows  that  the  problem  is  particularly  severe  near 
einc  =  11/1/2  (grazing).  On  the  other  hand  r(kz)  near  grazing 
is  well  known,  r(kz)  =  -1.  A  useful  way  to  look  at  the 
problem  is  to  represent  the  dispersion  relation  in  tne 
(kr  -  kz )  pi  ane  which  clearly  shows  the  increase  in  A kz  as  kr 
goes  from  0  to  kg. 
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(b),(c)  Sampling  in  the  z-domain 


The  sampling  interval  in  wavenumber  space,  Akz,  is  not 
chosen  through  aliasing  considerations  alone.  The  other 
consideration  is  the  resulting  resolution  in  the  z-domain 
imposed  by  the  total  sampling  interval  T, 

2n 

Az  =  — 

where  T  =  N  Akz 

A  decrease  in  Az  can  be  effected  either  by  a  decrease  in 
the  sampling  interval  in  the  wavenumber  domain  Akz  or  by  an 
increase  in  the  number  of  points  N.  It  is  easier  to  resort  to 
the  latter  method  as  the  reflection  coefficient  due  to  its 
rapid  decrease  for  large  vertical  wavenumbers  can  be 
conveniently  padded  with  zeros.  Moreover,  the  size  of  Akz, 
which  can  be  varied  in  a  synthetic  experiment,  is  usually 
fixed  in  an  actual  experiment. 

An  examination  of  the  reconstructions  of  a  half-space  of 
profile,  Figs.  14-17,  reveals  that  excellent  reconstructions 
can  be  achieved  for  appropriate  choices  of  N  and  Az.  A 
window  reconstruction  is  shown  in  Fig.  24.  The  main  effect 
of  a  decrease  in  az  is  that  adequate  reconstruction  of  the 
acoustic  profile  is  possible  to  a  greater  depth.  It  is  also 
true  that  for  given  sampling  intervals  Akz  and  az,  an 
increase  in  the  height  of  the  velocity  jump  at  the  seafloor 
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results  in  a  deterioration  of  the  reconstruction  with 
depth.  That  is  evident  for  a  half-space  by  comparing  Fig. 
14b  and  Fig.  17  and  for  the  Improved  Born  approximation  by 
analyzing  eg.  (6.18). 


6 . 4  Frequency  Scaling  of  the  Reflection  Coefficient 

Although  the  experiment  we  are  analyzing  is  mono¬ 
chromatic/  it  is  important,  for  a  proper  choice  of  operating 
frequency,  to  study  the  behavior  of  the  plane-wave  reflection 
coefficient  with  frequency.  For  the  simplest  acoustic 
medium,  a  homogeneous  half-space,  the  reflection  coefficient 
is  independent  of  frequency.  As  soon  as  a  spatial  scale  is 
introduced  in  the  acoustic  medium,  by  inserting  a  layer  for 
instance,  the  reflection  coefficient  becomes  frequency 
dependent.  This  dependence,  which  can  be  expressed  through 
the  continued  fraction  expression  of  eq.  (3.25),  is  at  the 
heart  of  the  Ware  and  Aki  inverse  method.  In  this  section  we 
show  that  the  solution  of  a  high  frequency  problem  is 
equivalent  to  the  solution  of  a  scaled  problem  at  a  lower 
frequency. 

The  one-dimensional  SchrOdinger  equation, 


(d  +  k2^(z)  =  V(z)<t>(z) 
dz  z 


(6.31) 


gives  rise  to  the  reflection  coefficient 


/ 


.  ik  z 
e  z 


2lk 

z 


V(z  '  )$  (z  '  )dz  ' 


(6.32) 
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At  a  different  angular  frequency  u>  '  =  atug,  with 
corresponding  water  wavenumber  k'  =  dkg*  the  Schr'Odinger 
equation  becomes 


(—  2  +  a 2k  2 )  <}>  ’  (z  )  =  V1  (z)$  '(z) 


(6.33) 


where  V'  represents  the  potential  at  the  new  frequency  u  ' . 
Now,  the  potential  is  frequency  dependent, 


Vz>  =  ko( 1  -  "2) 


(6.34) 


in  which  n,  the  index  of  refraction  is  a  function  of  depth  z. 
Therefore, 


V'  (z)  =  a2V0(z) 


(6.35) 


yields  the  Schrodinger  equation 


2 

+  a2kz2)^'  (z)  =  a2V0(z  )*' (z)  . 


(6.36) 


The  change  of  variable  z'  =  az  restores  the  original 
Schrodinger  equation  (6.31)  albeit  with  a  "stretched"  version 
of  the  original  potential, 


+  k 2 )  $  '  ( z  '  )  =  VQ  (z  1  /a  )<fr  '  (z  '  )  . 


(6.37) 
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The  corresponding  reflection  coefficient  is  now 


R'(kz)  =  /  TIF-  vo(zVaU'(z,)dz'  (6.38) 

-00  2 

The  reflection  coefficients  at  the  experimental 
frequencies  and  u '  are  therefore  identical  for  a  given 
wavenumber  kz  if 

VQ(z)  =  V0(z'/o)  •  (6.39) 

One  can  therefore  conclude  that  a  high  frequency  experi¬ 
ment  (u*  >  u)  is  equivalent  to  a  low  frequency  experiment 
with  a  stretched  profile  (a  >1).  One  can  therefore  compare 
the  reconstruction  of  a  given  profile  at  two  frequencies,  say 
a)  and  oj '  =aa i,  by  comparing  the  reconstruction  at  a  single 
frequency  a  of  the  given  profile  with  its  stretched  version 
(stretch  factor  a).  Now,  it  is  a  numerical  fact,  as  seen  in 
previous  examples,  that  profile  reconstruction  via  the 
Gelf and-Levitan  algorithm  deteriorates  with  depth.  It  is 
therefore  clear,  at  least  for  simple  profiles,  such  as  a  step 
(half-space)  or  a  window  (layer),  that  a  lower  experimental 
frequency  entails  deeper  reconstruction  as  shown  in  Fig.  30. 
That  holds  for  more  complicated  profiles,  and  we  therefore 
conclude  from  a  frequency  scaling  point  of  view  that  the 
lowering  of  the  experimental  frequency  allows  for  deeper 
profile  reconstruction  without  noticeable  effects  on  the 
reconstruction  of  the  detailed  variations  of  the  profile. 
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6 . 5  Profile  Reconstruction  in  the  Presence  of  Density 


Variations 

Although  the  focus  of  this  thesis  has  been  on  the 
reconstruction  of  the  acoustic  profile  in  the  presence  of 
velocity  variations  at  a  single  frequency,  it  is  useful  here 
to  show  how  one  can  extend  the  method  to  the  retrieval  of  the 
velocity  in  the  presence  of  density  variations  and  also  to 
recover  the  density  profile. 

We  have  shown  in  Chapter  II  (eq.  (2.5))  that  in  the 
presence  of  smooth  density  variations,  a  change  of  variables 
retained  the  governing  Helmholtz  and  SchrOdinger  equations 
with  an  attending  redefinition  of  the  index  of  refraction  to 
account  for  variations  in  the  density  p , 


=  n2  +  k 


-2 


v2p 


3r  1  „  ,2, 

~  Tlr  7p)  ) 


(6.40) 


The  density  dependent  potential  to  be  reconstructed  is  now, 


(6.41) 


A  single  frequency 
n ' (z ) .  To  recover 
and  p(z)  one  needs 
frequencies  and 
respectively.  The 


experiment  can  only  hope  to  reconstruct 
n(z)  (and  therefore  the  velocity,  c(z)) 
to  carry  out  the  experiment  at  two 
u2  with  water  wavenumbers  kg  and  kj, 
associated  potentials  are  then 
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V t ( z )  =  k\{l  -  n'J) 


and 


(6.42) 


V2(z) 


These  potentials  can  be  reconstructed  by  the  application  of 
the  Gelfand-Levitan  algorithm  to  the  corresponding  reflection 
coefficients.  One  can  then  obtain  the  difference  of  the 
potentials , 


Vl" 


V9  "  (k1 


(6.43) 


and  therefore  retrieve  the  velocity  dependent  index  of 
refraction  n(z). 


n(z ) 


1 


(6.44) 


and  subsequently  reconstruct  the  acoustic  velocity  profile, 


c  (z ) 


_  C0 

'  nTzT  * 


(6.45) 


One  of  our  numerical  computations,  was  to  conduct  the 
two-frequency  procedure  on  the  profile  of  Figure  10.  The 
reconstruction  of  c(z)  is  shown  in  Figure  3 0 i .  The  recon¬ 
struction  features  in  this  case  of  c(z)  are  similar  to  the 
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constant  density  case.  The  density  as  function  of  z  could  in 
turn  be  recovered  by  solvinq  the  differential  equation  (6.40) 
for  p(z)  qiven  n(z). 

Density  discontinuities  invalidate  the  asymptotic 
behavior  of  the  reflection  coefficient  presented  in  Chapter 
III.  In  fact,  the  reflection  coefficient  is  generally  not 
integrable,  and  one  has  to  introduce  generalized  functions  in 
its  Fourier  transform.  The  attending  difficulties  and  their 
resolution  in  the  Ge If and-Levitan  algorithm  have  not  been 
studied  here. 

6 . 6  Acoustic  Attenuation 

The  study  of  the  attenuation  of  acoustic  waves  in  marine 
sediments  has  been  studied  recently  by  Rajan  and  Frisk^67^ 
who  proposed  a  perturbative  inverse  method  for  the  recovery 
of  the  attenuation  data  from  reflection  data.  Rajan  has  suc¬ 
cessfully  inverted  for  the  acoustic  attenuation  profile  given 
the  reconstructed  velocity  profile  we  had  obtained  through 
the  Gelf and-Levitan  algorithm ^ .  Here,  we  shall  look  at 
the  effects  of  intrinsic  attenuation  on  the  reconstructed 
acoustic  profile. 

Although  our  formulation  of  the  model  of  a  Frisk  experi¬ 
ment  does  not  include  attenuation,  it  is  possible  to  posit  a 
lossy  acoustic  profile,  generate  the  corresponding  reflection 
coefficient  and  then  run  it  through  the  Gelfand-Levitan 


-132- 


algorithm.  As  expected,  for  small  loss  (a  ~  0.001-0.005 
dB/m),  the  reflection  coefficient,  its  Fourier  transform,  and 
the  reconstructed  profile  are  little  affected  by  the  per¬ 
turbation  (see  Figs.  25a,  25b,  26,  27).  As  loss  increases 
(a  ~  0.01  dB/m),  the  reconstructed  profile  deteriorates 
rapidly.  It  should  be  pointed  out  that  the  lower  values  of 
intrinsic  attenuation  prevail  in  the  sediments  in  the  abyssal 
plain  at  the  frequencies  of  interest  (220  Hz).  Incidentally, 
one  of  the  advantages  of  a  monochromatic  experiment  is  that 
the  frequency  dispersion  law  of  the  intrinsic  attenuation, 
which  is  difficult  to  establish  experimentally,  particularly 
at  low  frequencies  (<  1  kHz),  does  not  enter  into 
consideration. 

6 . 7  Noise 

As  demonstrated  in  Chapter  V,  the  Gelf and-Levitan 
algorithm  is  stable,  with  small  errors  in  the  reflection 
coefficient  resulting  in  small  errors  in  the  reconstructed 
velocity.  The  numerical  experiments  we  have  conducted  by 
adding  zero-mean  Gaussian  stationary  noise  to  the  reflection 
coefficient  input  support  our  previous  conclusion. 

Gaussian  noise  was  added  to  both  the  real  and  imaginary 
part  of  the  reflection  coefficient  generated  by  the  method  of 
Section  6.1.  The  resulting  reflection  coefficient  was 
submitted  to  the  usual  steps  involved  in  the  inversion 
procedure,  namely  the  Fourier  transform  step  followed  by  the 
application  of  the  Gelfand-Levitan  algorithm. 


-133- 


The  signal  to  noise  ratio  (SNR)  is  defined  as 


SNR  =  10  log^ty^-) 


where  the  signal  and  noise  powers  for  N  discrete  points  are 


1  N  2 
rs  -ff  I,  K' 


n  =  l 


and 


2 

I  =  a  i  the  variance  of  the  noise, 
n  n 

As  shown  in  Fig.  31  for  a  window  profile,  the  perturbat¬ 
ion  of  the  reflection  coefficient  by  the  addition  of  zero- 
mean  Gaussian  noise  (a  =  0.1)  results  in  a  roughly  propor¬ 
tional  degradation  of  the  reconstructed  potential. 

The  preliminary  assessment  of  the  effect  of  noise  leads 
us  to  conclude  that  the  Gelfand-Levitan  inverse  method  is 
stable  in  the  presence  of  noise.  This  analysis  can  be 
refined  in  the  future  by  including  a  more  pertinent  model  for 
the  noise  based  on  Mook '  s  results^^;  Mook  has  shown  that 
the  addition  of  zero-mean  stationary  white  Gaussian  noise  to 
the  point  source  pressure  field  resulted  in  a  non-stat ionary 
variance  with  the  noise  power  concentrated  near  the  vertical 
(kz  =  kg,  the  water  wavenumber).  One  should  also  note  that 
the  measured  reflection  coefficient  itself  can  be  regarded 
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as  a  non-stationary  process  with  a  mean  that  varies  between 
one  near  grazing  incidence  and  zero  at  infinity.  That 
consideration  may  lead  to  a  more  appropriate  measure  of 
performance  than  that  of  SNR  presented  here. 
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CHAPTER  VII 


CONCLUSION 

In  this  thesis  we  have  examined  the  theoretical  and 
computational  underpinnings  of  a  novel  approach  to  the 
determination  of  the  acoustic  parameters  of  the  oceanic 
sediment  layer.  Traditional  marine  seismic  methods,  acoustic 
reflection  and  refraction  measurements,  "yield  no  velocity 
information  in  the  top  of  the  first  sediment  layer  which  is 
of  critical  interest  for  modeling  the  sea  floor  for  under¬ 
water  acoustics"  (34  ).  That  is  precisely  where  the  method 
we  have  analyzed  in  the  thesis  is  most  accurate.  Interval 
velocity  calculations  are  restricted  to  layer  thicknesses 
larger  than  one  twelfth  the  water  depth  (34  ),  and  subsequently 
the  velocity  at  the  sediment  interface  cannot  be  determined 
accurately.  In  the  absence  of  in  situ  core  measurements 
of  surface  velocity  VQ ,  the  velocity  gradient  at  the  top 
of  the  sediment  column  is  also  uncertain.  The  direct 
inverse  method  presented  here  requires  only  the  a  priori 
knowledge  of  the  speed  of  sound  in  the  water  above  the 
sediment  layer. 

By  operating  at  a  single  frequency,  the  effects  of 
dispersion  are  separated  from  the  propagation  process. 

Moreover,  the  dispersive  characteristics  of  the  medium  can 
be  studied  by  performing  the  experiment  at  various  frequencies. 
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A  single  frequency,  steady-state  experiment  also  means 
that  the  measurement  relies  on  amplitudes  rather  than  on 
arrival  times  (explosive  methods).  While  it  is  true  that 
time  is  measured  more  accurately  than  acoustic  pressure, 
arrival  times  can  not  always  be  interpreted  correctly  due 
to  multiple  reflections.  Amplitude,  in  turn,  may  be  affected 
by  a  host  of  factors  to  which  arrival  time  is  insensitive 
such  as  loss  and  diffraction;  it  is  not  yet  possible  to  judg~ 
the  relative  merits  of  the  two  methods  in  the  absence  of  appro¬ 
priate  experimental  data. 

Contributions 

The  determination  of  the  acoustic  properties  of  the 
ocean  bottom  was  shown  to  be  equivalent  to  the  reconstruction 
of  an  unknown  potential  in  a  Schrodinger  equation  from  the 
plane-wave  reflection  coefficient  given  at  all  angles  of 
incidence  (Chapter  II) . 

The  pivotal  role  of  the  reflection  coefficient  lead 

us  to  a  detailed  examination  of  its  properties  in  Chapter  III. 

In  particular,  we  showed,  by  induction,  that  the  reflection 

2 

coefficient  decays  at  least  as  rapidly  as  (1/k  )  and  is 

therefore  integrable. 

The  derivation  of  our  approach  to  the  direct  inversion 
method  was  presented  in  Chapter  IV.  The  Gelf and-Levitan 
method  was  extended  to  the  case  where  the  acoustic  velocities 
on  either  side  of  a  slab  are  different.  That  is,  of  course, 
the  case  in  the  ocean  bottom  problem  where  the  acoustic  velocity 
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in  the  basement  is  larger  than  the  acoustic  velocity  in  the 
water.  We  also  showed  that  the  neglect  of  bound  states  is 
justified  at  the  current  operating  frequency  in  both  clayey 
silt  and  in  silty  clay. 

Three  methods  for  the  numerical  solution  of  the 
Gelfand-Levitan  integral  equation  were  investigated  (Chapter 
V) .  The  first  method  we  developed  is  a  series  expansion  of 
the  solution  obtained  by  successive  approximations.  The 
first  two  terms  of  this  expansion  represent  a  substantial 
improvement  over  the  well  known  Born  approximation. 

The  other  two  numerical  methods  presented  in  Chapter  V 
are  based  on  the  discretization  of  the  Gelfand-Levitan 
integral  equation.  They  represent  two  ways  to  bypass  the 
matrix  inversion  inherent  in  a  straightforward  solution  of 
the  discretized  equation.  We  then  obtained  estimates  for 
the  bound  on  the  error  in  the  integral  of  the  potential  due 
to  discretization  errors  and  due  to  errors  in  the  reflection 
coefficient. 

In  Chapter  VI,  we  discussed  the  numerical  results 
obtained  from  the  inversion  of  synthetic  data.  By  dealing 
with  synthetic  data,  we  insured  that  the  bottom  profile  was 
known  exactly  and  that  the  effectiveness  of  the  method  could 
be  studied  without  any  fear  of  experimental  imperfections. 

We  concluded  that  the  Gelfand-Levitan  method  appears  to  be 
very  accurate  at  the  top  of  a  sediment  column,  just  where 
the  more  usual  methods  are  least  accurate.  The  resolution 
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obtained  is  less  than  the  wavelength  of  the  acoustic  source 
in  the  water.  The  degradation  of  the  reconstructed  velocity 
profile  becomes,  however,  pronounced  if  the  reflection  data 
is  restricted  to  real  angles  above  critical.  Perturbations 
of  data  were  also  studied.  Perturbations  such  as  intrinsic 
loss  in  the  acoustic  medium  cr  noise  in  the  data  produce 
perturbations  in  the  reconstructed  profile.  The  inclusion 
of  density  variations  requires  the  use  of  two  frequencies 
and  two  separate  inversions  of  the  Gelf and-Levitan  equation. 

We  were  able  to  gauge  the  performance  of  the 
numerical  schemes  through  the  study  of  the  inversion  of  the 
acoustic  profile  for  two  half-spaces  (constant  velocity) , 
for  which  we  derived  the  Fourier  transform  of  the  reflection 
coefficient  analytically.  The  improvements  wrought  by 
the  improved  Born  method  are  clear,  as  are  the  effects  of 
sampling  on  the  reconstructed  profile.  The  improved  Born 
method  represents  a  fast  and  easy  to  implement  method  of 
recovering  the  velocity  at  the  top  of  the  sediment  column. 
The  two  finite  difference  methods  are  more  time-consuming 
but  yield  an  accurate  reconstruction  of  the  acoustic  profile 
over  greater  depth. 

Future  Work 

On  the  theoretical  front,  we  suggest  a  thorough 
investigation  of  the  Gelf and-Levitan  method  in  the  presence 


139 


of  density  variations.  It  may  also  be  useful  to  incorporate 
loss  directly  into  the  original  formulation.  One  should  also 
seek  efficient  numerical  implementations  of  the  Gelfand- 
Levitan  algorithm  that  could  increase  the  penetration  depth 
of  the  reconstruction.  In  this  respect,  we  think  that  the 
combination  of  the  Gelf and-Levitan  algorithm  with  a  priori 
information  such  as  the  acoustic  properties  of  the  basement 
might  constitute  a  promising  approach. 

On  the  experimental  front,  the  testing  of  the  Gelf and- 
Levitan  inverse  method  on  field  data  should  be  given  priority 
to  determine  its  ultimate  value.  If  an  actual  ocean-based 
experiment  were  precluded  at  the  moment,  we  would  suggest 
carrying  out  a  similar  electromagnetic  experiment  on  dielectrics 
at  microwave  frequencies. 

On  the  numerical  front,  one  would  want  to  test  the 
whole  experimental  scheme,  starting  from  pressure  measurements 
due  to  a  point  source  and  ending  with  the  Gelfand-Levitan 
inversion.  Special  attention  should  be  paid  to  the  effect 
of  noise  in  the  recording  of  pressure  on  the  plane-wave 
reflection  coefficient  and  ultimately  on  the  reconstructed 
acoustic  parameters. 
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It  is  remarkable  that  the  Frisk  experiment  which  was 
conceived  on  purely  intuitive  grounds,  when  modelled 
theoretically,  bears  out  the  expectation  that  accurate  results 
are  achievable.  To  date,  experimental  data  has  been  inter¬ 
preted  by  time  -  consuming  trial  and  error  procedures.  Our 
mathematical  and  numerical  approach  suggests  chat  a  more  direct 
inverse  method  for  processing  experimental  data,  requiring 
no  a  priori  information  about  the  acoustic  parameters  of 
the  bottom  is  feasible. 
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Figure  11:  Frisk.  Profile 
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VERTICAL  UAVENUMBER 

Figure  12a:  Reflection  Coefficient  Amplitude  for  a  Half-Space 
f  =  25Hz 
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DEPTH  IN  METERS 

Figure  14a:  Reconstruction  for  a  Half-Space  (Potential) 
Az  =  0.767m,  f  =  25Hz 
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Figure  14b:  Velocity  Reconstruction  for  a  Half-Space  (m/s) 
R(z)  obtained  numerically,  iz  =  0.767m 
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Figure  17a:  Velocity  Reconstruction  for  a  Half-Space. 

cn  =  1500m/s ,  c.  =  1600m/s) .  Az  =  0.3125m 
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Figure  18c:  R(z) ,  Fourier  Transform  of  Figure  18 (a, b) 
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Figure  18d:  Frisk  Profile  Reconstruction  -  Nq  Iteration 
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Figure  18e:  Frisk  Profile  Reconstruction  -  20  Iterations 
Az  =  0.767m,  f  =  25Hz 


Frisk  Profile  Reconstruction  -  40  Iterations 
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DEPTH  IN  METERS 

Figure  18g:  Velocity  Reconstruction  -  Frisk  Profile  -  (Detail) 
Az  =  0.767m,  f  =  25Hz 
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Figure  19d:  Reconstructed  Velocity  -  Frisk  Profile 
Limited  Range  K  =  0.512,  f  =  25Hz 
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Figure  20:  Frisk  Profile  Reconstruction.  50  Iterations 
AK  =  0.0005.  Az  =  3.06m,  K  =  0.125 
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DEPTH  IN  METERS 

Figure  23:  Frisk  Profile  Reconstruction. 

Limited  Range  K  =  0.016,  A z  =  3.06m 
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Figure  26a:  Reflection  Coefficient  Amplitude.  f  =  25Hz 
Lossy  Half  Space,  ct  =  O.OOOldB/m 
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VERTICAL  UAVENUMBER 

Figure  27a:  Reflection  Coefficient  Amplitude, 
Lossy  Frisk  Profile,  a=  0.005dB/m 


VERTICAL  UAVENUMBER 

Figure  27b:  Reflection  Coefficient  Phase,  f  =  25 
Lossy  Frisk  Profile,  a=  0.005dB/m 
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Figure  27d:  Velocity  Reconstruction ,  f  =  25Hz 
Lossy  Frisk  Profile,  a-  0.005dB/m 
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Figure  28:  Stickler  Profile  Reconstruction.  Az  -  6.03m,  L  =  40m,  f  = 
1:  Original  Profile,  2:  Finite  Difference,  3:  Improved  Born 
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Figure  29:  Stickler  Profile  Reconstruction.  Az  =  6m,  L  =  90.6m 

1:  oroginal  profile,  2:  Finite  Difference,  3:  Improved  Born 
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Figure  30i:  Reconstructed  Velocity  Potential  V(z)  -  p (z)  =  c(z)/1000 
1:  Original  Stickler  Profile,  3:  Improved  Born 
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